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ABSTRACT 
INVESTIGATIONS OF PRIMAL TECHNIQUES FOR INTEGER LINEAR PROGRAMS 


by 
MARK ANDREW HOPPER 


Submitted to the Department of Civil Engineering on 17 May 1968 in 
partial fulfillment of the requirements for the degree of Civil Engineer. 


A primal feasible (all-integer) integer linear programming algo- 
rithm, developed independently by R. H. Gonzales and R. D. Young, has 
been coded for a digital computer and has been implemented as a problem- 
oriented language under the ICES/OPTECH subsystem for mathematical pro- 
gramming. Complete program documentation and a guide for the use of 
this system are included. 


Some convergence problems of the Gonzales-Young algorithm are 
identified and geometric interpretations that help to illuminate these 
problems are offered. 


A strong relationship between the Gonzales-Young algorithm and 
R. E. Gomory's Algorithm I (a dual feasible method) is established. 
This relationship and the interpretation of convergence problems in 
the Gonzales-Young algorithm are combined to produce a modified, and 
hopefully more efficient, primal algorithm. Performance of this algo- 
rithm is measured against other integer linear programming algorithms 
on a small set of test problems. 


The concepts of the modified primal algorithm for the all-integer 
problem are enlarged upon to develop the essentials of the first primal 


mixed integer-continuous linear programming algorithm that has ever been 
proposed. No computational results of this algorithm are available. 


Thesis Advisor: Dr. Alan ?'. Hershdorfer 


Title: Assistant Professor of Civil Engineering 
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CHAPTER | 


INTRODUCTORY REMARKS 


1.1] Integer Linear Programming Problems 


The formulation of an Integer Linear Programming (ILP) problem 
is similar to that of the more widely known Linear Programming (LP) 
problem. Like the LP formulation the ILP formulation is character- 
ized by a linear objective function that is to be optimized, subject 
to a set of linear convexity constraints and non=-negativity restrice- 
tions on the problem variables. However, the ILP formulation imposes 
the additional restriction that some or all of the problem variables 
must take on integer values in a feasible solution. 

The ILP formulation arises quite frequently in both business and 
engineering, Integrality restrictions on the problem variables some= 
times result from the inherent indivisibility of the units that these 
variables represent, such as men and machines, At other times inte=- 
grality is employed in the formulation in a somewhat artificial manner 
Since it can be used conveniently to express logical relationships 
(e.g. decision variables), Regardless of how the integrality restric- 
tions arise, special procedures are necessary to ensure that they are 


satisfied. 





Since the ILP formulation is so similar to the LP formulation, it 





was a natural consequence that the first attempts to solve the ILP 
problem centered around the already proven Simplex method for linear 
programming. The first such approach was forwarded by G. B. Dantzig(1). 
This technique attacked the All Integer Linear Programming (AILP) prob- 
lem by successively modifying the linear program that was inherent in 
the AILP formulation (if the integrality restrictions were ignored). 
The key aspect of Dantzig’s approach was that he deduced new, more 
binding, constraints from the optimal tableau of the linear program. 
These additional constraints did not exclude any feasible integer 
points, but they made the available non-integer optimal extreme point 
infeasible, Thus the addition of such constraints to the previous 
constraint set resulted in a reduction of the solution space and re= 
optimization was required, The supposition was that by successively 
reducing the solution space, while not excluding integer solutions, 
eventually one of the reoptimized tableaus in the sequence would have 
to display the desired integrality properties in its solution (2] . 
Although counterexamples ultimately showed that Dantzig’s technique 
could not guarantee a finite convergence, the concept helped to launch 
the very important class of solution techniques for ILP problems that 
are collectively referred to as the “cutting plane" methods. The 
commonality among these techniques arises from the fact that they al] 
employ additional constraints (cuts) that function in a fashion similar 
to those proposed by Dantzig’s approach. 

The first proven algorithms for solving ILP problems by the cutting 
plane approach were developed by R. E. Gomory [3,4,5}. The cuts that 


A 





are generated in these algorithms are founded on a much more 
sophisticated mathematical basis than the proposals of Dantzig. In 
fact, while Dantzig can be credited with proposing the first legiti- 
mate cut, Gomory’s work continues to provide the basis for most of 

the meaningful developments within the cutting plane context. Gomory’s 
first algorithms have been followed by numerous investigations into 

the nature of the various cuts that can be generated in his algorithms. 
The manner in which cuts are chosen from the many that are available 
at a given stage of the algorithms has been shown to affect the per- 


formance of the various algorithms to a considerable degree |6]. 


1.3 The Motivation for a Primal Algorithm 





Even today there is no single algorithm that can be said to handle 
all ILP problems in a satisfactory manner. This is evidenced by the 
fact that the topic continues to enjoy a good deal of attention in 
technical publications. Gomory’s algorithms and the various modifica- 
tions to them have enjoyed only sporadic success, Some problems have 
been solved in what has been considered a reasonable amount of time 
while others have proved computationally expensive beyond tolerable 
limits, Unhappily the computation required for the solution process 
has not always been found to be proportional to the size of the problem, 
so that no reliable guidelines are available for predicting the cost of 
solving a particular problem, This difficulty has been compounded by 


the fact that, until recently, all of the cutting plane algorithms have 
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been developed as “dual feasible” methods. This characteristic meant 
that no primal feasible solution to the ILP was available until the 
optimal solution was found, As a consequence all of the problems that 
proved too costly to solve to optimality resulted in a total waste of 
computational effort. Thus the need for a "primal feasible" method was 
evident. Such a method would provide at least a feasible solution to 
the problem under consideration. This solution, though not proven to 
be optimal, could be useful and sometimes might be fairly close to the 


optimal. 


1.4 Contributions to the Primal Algorithm 
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In June of 1965 Gonzales [7] and Young [8] reported independent 
development of the first primal feasible algorithms designed to solve 
the AILP problem. The two algorithms are basically the same in that 
they are both direct extensions of proposals made by Ben=Israe] and 
Charnes [9]. The only substantial difference between the two algorithms 
lies in the tableau representation and not in the computational approach. 
Both of these algorithms use a cut that was developed by Gomory in what 
is referred to as his Algorithm II [6). In fact, at least in their 
basic forms, these algorithms represent a modification of Gomory’s 
Algorithm II whereby primal pivot selection rules are substituted for 
dual pivot selection rules. This does not in any way detract from the 
importance of either Young’s or Gonzales’ work. Since Gomory did not 


address himself to the applicability of his cuts to a prima! approach 
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the credit for the primal algorithm belongs to Young and Gonzales. 
Young showed that his version of the algorithm could be made to con- 
verge in a finite number of iterations. Although Gonzales was not 

able to prove finite convergence save for the two dimensional problem, 
he did provide a clever method for achieving a starting basic feasible 
solution when one is not available in the initial tableau. In order 
to give dual credit for the development of the primal algorithm, the 
two versions will be referred to here as one and will be called the 


Gonzales-Young (G=Y) algorithm [10). 


1.5 Scope of the Current Study 


The work that is reported in this paper was carried on in two 
related parts, The primary objective of this research has been to 
provide for ICES/OPTECH Subsystem [I1]a useable capability for solving 
AILP problems in the form of the Gdnzales=Young algorithm. This capae 
bility will not only be valuable to researchers in this particular 
branch of mathematical programming, but will also provide a useful 
tool for research in civil engineering and management in general. 

At the same time that work was being carried forward in the come 
puter implementation phase, the mechanics of the G=-Y algorithm were 
being investigated with a view toward improving the convergence 
characteristics of the algorithm itself. Originally this second 
aspect of the research was intended to uncover an optimum pivot 
selection criteria and perhaps to discover ways of generating cuts that 
would be more binding than the ones that are generated in the G-Y 


algorithm. In this regard the author was able to find some primal cuts 
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that were better than the normal G-Y cuts. Unfortunately each of 

these new cuts resulted in a loss of the all-integer tableau. As a 
result these cuts could not be used in the G-Y context of all integer 
solutions. However, the demonstrated power of the different cuts led 
the author to attempt a modification of the G-Y algorithm. The motiva- 
tions behind this modification and the results that have been attained 
are reported in the second part of this paper. Finally, the basic ap- 
proach of this modified AILP algorithm is extended to produce a primal 
algorithm for the MILP problem. 

To summarize, the research reported in this paper has been pre- 
sented in two parts. The first part describes in detail the consider- 
ations that went into the computer implementation of the G-Y algorithm 
(called OPAILP). This part has been written in such a way that it can 
Stand alone as a user's guide for the OPTECH all integer linear pro- 
gramming capability. The second distinct part includes an interpre- 
tation of some of the convergence problems of the G-Y algorithm that 
served to motivate the proposal of a modified AILP algorithm. This 
discussion is extended to relate the G-Y algorithm to Gomory's Algorithm I 
for the AILP problem. A modified algorithm is proposed that is essen- 
tially a hybred of the G-Y algorithm and Algorithm I. Like the G-Y 
algorithm the modified algorithm is primal feasible, but it specifically 
avoids the difficulties that have been shown to impede the G-Y conver- 
gence process. The performance of this hybred algorithm is measured 
against several AILP algorithms (including G-Y) on several problems. 
Finally, an outline of the essentials of a primal MILP algorithm are 


given, though no computational results are reported. 
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CHAPTER II 
THE GONZALES-YOUNG AILP ALGORITHM 


In order to provide a reasonable basis for subsequent discussions 
of the computer implementation of the G-Y algorithm, this chapter will 
outline the fundamental steps that the algorithm performs. The prob- 
lem formulation and notation employed here were used by Gomory in his 
first ILP algorithm and by Gonzales in his presentation of the primal 
AILP algorithm. This approach was chosen over that of Young because it 
facilitates comparisons among the ILP algorithms. These comparisons are 


undertaken in the second part of this paper, mainly in Chapter IV. 


2.1 The Problem: Notation and Formulation 
The G-Y primal AILP algorithm will solve a problem that can be 


formulated in the following way: 


Find non-negative integers Xs 576 oil 


that maximize 


| 
» CX, (2.1.1) 


Subject to 


Die ade eth (2a) 


as 
a5 4X4 <b, 


Cefe 


iJ 


cs. 
ii 
om) 


with the Cas diye and D. integers 


Sas 





Except for the integrality restrictions on the problem variables this 
formulation is exactly that of the classical LP formulation. This 
problem notation is altered slightly by Gonzales in order to recast 
the problem into a parametric form suggested by A. W. Tucker.. 

Each of the inequality constraints (2.1.2) is made into an equation by 
adding a slack variable Si. Since all of the constraint coefficients 
are integer and the Xs are restricted to non-negative integers, the Ss. 
are also restricted to non-negative integer values for a feasible 
solution. Non-negatively constrained parametric variables t are 
introduced to act as proxies for the original problem variables. For 
notational consistency the Cc; are referred to as a9 j and the b. as 
aio: Thus in terms of the Tucker notation the problem can be restated 


in the following form: 


maximize 2 


n 
when: Z=ay + ) “ap j(-t,) (2s) 
j=l 
n 
s;7a; 9 + » a; s(-t,) fel... (2.1.4) 
j=l 
x. (-1)(-t,) j=l en (Zalie5) 


where the te are restricted to non-negative integer values and z is an 


unsigned variable. 


ee 





The above relations can be expressed in a tableau format as 


follows: 





Figure 2.1.1 


In this tableau the problem variables are listed along the left margin, 
The solution that corresponds to any given tableau is found in the 
zeroth column. The problem variables may be viewed as slack variables 
associated with the hyperplane representation of the problem constraints 


in non-basic (parametric) space. 


ae 





2.2 Pivoting in the Tableau 


By employing Gaussian elimination the problem variables can be 
expressed in terms of a new set of non-basic (parametric) variables. 
Completion of this process is referred to as a pivot step. The com- 
putations that are required in a pivot step are referred to as a 
tableau update. Pivot selection consists of finding a parametric 
variabl ta that can be raised to a positive level to increase the 
value of the objective function and finding the constraint row that 
most strictly binds the increase of t,- From examination of the 
expression (2.1.3) it is evident that the objective value can in- 
crease with an increase in t only if 49.4 < 0. This establishes a 
criteria for selecting a possible pivot column. Any column j for which 
99 3 < 0 may be selected for pivoting. Likewise a problem variable Sy) 
can decrease toward infeasible values with an increase in t only if 
ak a > 0. Thus for any column gq the most binding constraint is found 
by: 

“a (a; 9/8; 4) i=1,...,mtn ne 0 (252e) 
(a/b denotes a divided by b) 


Having selected a pivot column q and a pivot row k, the tableau element 


ay q is referred to as the pivot element. The following formulae update 
the tableau: 
column q Sena “a; /ay g) 1=0,1 5.2 0gntth (2.2.2) 
columns jfq 95 5-94, 57 94 962k, 5/2k a? i=0,1,...,mtn (2.23) 
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Geometrically the parametric variables may be viewed as an 
n-dimensional orthogonal coordinate system. Since the parametric 
variables are restricted to take on only non-negative values, the 
optimal solution to the problem always lies in the positive orthant 
of this coordinate system. The updating operations translate the 
origin of the parametric coordinate system from one vertex of the con- 
vex polyhedron to an adjacent one. The coordinate axes also undergo 
a rotation such that the convex polyhedron is continually distorted 
as a result of the pivot operations. The current solution is the 
vector of values found in the zeroth column when the parametric vari- 
ables are at zero levels. Hence the current solution always lies at 
the origin of the parametric coordinate system. At optimality a co- 
Ordinate system with its origin at the optimal extreme point of the 


convex polyhedron has been attained. 


2.3 The Gomory Cuts 

It is evident that the updating formulae (2.2.2, 2.2.3) could not 
always be expected to result in integer solutions to the problem. 
However, by examining these formulae one can deduce that such will be 
the case if all of the entries in the initial tableau are integer and 
if all pivot elements (ay ig in 2.2.2 and 2.2.3) are selected as ‘11. 


These requirements are met by the special Gomory cuts that are employed 


in the G-Y algorithm. 


sje 





In terms of the tableau representation, Gomory consicered a 


constraint of the type: 


n 
aes +) Ben gkstit (2.3.1) 
a 
and showed that non-negative integer vectors t = (t.....st,) which 


Satisfy (2.3.1) also satisfy: 


n 


r= = [a, o/pl +) La; ;/p](-t,) (253.2) 
j=l 


p 0 


and [a/p] denotes the largest integer less than or equal to a/p. 


The G-Y algorithm considers the most binding constraint k for a ‘column q 
to be (2.3.1). It then selects P=ay (the pivot element) so that the 

cut that is generated will have +1 in the pivot column. The cut (2.3.2) 
is appended to the tableau in row mtn+l and the pivot element is changed 


from a +]. Since the initial tableau has only integer 


k,q °° 4nentl.q ~ 
entries, and since the coefficients of (2.3.2) are all integers, in- 
tegrality is maintained in the updating process. The following is an 
example of one pivot step on the problem: 


max Zz = 2x4 + 3x 
5 
5 
X12 Xo non-negative integers 


-Xy + exo 


IA FEA fo 


on ee 


ale 





The initial tableau is shown in Figure 2.3.1]. The pivot element of 
Figure 2.3.1] was a5 3 = 2°. The cut r was generated from constraint Sy. 


The new pivot element ae 3 = 1* produces an updated tableau shown in 


magure 2.3.2. 





Figure 2.3.1 Figure 2.3.2 


These operations are shown geometrically in Figures 2.3.3, 2.3.4, and 
2.3.5. Figure 2.3.3 shows the original problem constraints. Figure 2.3.4 
shows the cut that is added and indicates the move from integer vertex ] 
to integer vertex 2. Figure 2.3.5 shows the convex polyhedron in terms 

of the new parametric coordinate system (notice the distortion of this 


polyhedron). 


aCe 








Figure 2.3.3 
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Figure 2.3.4 





2.4 Some Properties of the Cuts 


If in the process of generating a cut (2.3.2) it occurs that 
Lay 9/ ay od = 0, the updating operations will not change the previous 
solution. Such a cut is referred to as a "zero cut" and the updating 
Operations that are involved in pivoting on such a cut are referred to 
as a "Stationary iteration". If Lay o/ ay gd # 0 then "breakthrough" 
is said to have occurred and the resulting update is referred to as a 
“transition iteration”. 

The fact that zero cuts can (and often do) arise in the course of 
the algorithm presents a problem of convergence not unlike that en- 
countered in the normal simplex method for linear programming. (Dis- 
cussion of this relationship will be postponed until Chapter IV). Young 
designed a method that takes the possibility of zero cuts into account 
and guarantees that the algorithm will converge in a finite number of 
pivot steps. Although Young's proof is too elaborate to reproduce here, 
it 1s apparent that zero cuts should be avoided wherever possible since 


they slow down the solution process. 


2.5 Summary of the Basic Gonzales-Young Algorithm 


All of the essentials of a basic G-Y algorithm have been presented. 
Although the convergence proof of the algorithm depends on some rather 
Specialized pivot selection rules, the basic algorithm can be presented 
in a straightforward manner. To date no problem has been shown to cycle 


under the basic algorithm if the problem can be guaranteed to be bounded 


= 





in every direction (i.e., possesses a closed convex set). The algorithm 


assumes an all integer tableau. 


A. Among the columns j with Aq 24 O select one for pivot 
consideration. (If there are no ay 3s 0, the current 


tableau is optimal.) 


B. Having chosen column q for pivot consideration, select 
the most binding row k from: 
OF 1=1 52... 


a 
(If there is no aig? O i=1l,...,mtn, the problem is un- 
9 


bounded. ) 
C. From row k generate a cut by setting p=a, : inmlaeo.e) 


p>” If La, 0/4 q = 0, check to see if there are other 
negative columns that have not been considered. If there 


are, go back to A. If not, go to E. 


E. Append the cut in row mtnt+l of the tableau and select 
antnt] a = 1 as pivot element. Perform the updating 


operations, erase the cut, and go back to A. 
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2.6 Phase I of the G-Y¥ Algorithm 


It sometimes occurs that the AILP problem cannot be formulated in 
such a way as to provide a starting basic feasible solution in the 
initial tableau. There are two forms of infeasibility that can arise in 
the initial tableau. The first case is when equations instead of in- 
equalities appear in the problem formulation. The second case is when 
the < inequalities have a negative right hand side (i.e. the constraint 
Slack has a negative value in the initial tableau) 

Suppose that constraint k is originally formulated as an equation. 
The slack Si. is still added to this constraint when it is changed to 
the tableau format. However, any feasible solution to the equation 
requires that S,=0. Gonzales shows that, by successively pivoting on 
cuts generated from the row marked Sp the row can be reduced to the 


representation s,= O-a mca if there is a feasible integer solution 


k 5q 
to the equation. At this point it is apparent that t, can be kept to 
a zero level simply by forcing it to remain non-basic. Therefore, the 
column marked “t, can be dropped from the tableau and subsequent pivot 
Steps assure that Sy. remains at a zero level (feasible). These oper- 
ations are applied to the equations one by one until all have been 
Satisfied. 

When any tableau has a negative problem variable (including the 
Slacks), this variable is currently infeasible. Gonzales reasons that 
the solution space associated with any subset of the problem constraints 


must include the solution space associated with the complete constraint 


set. Since the goal is to raise a currently negative slack variable Sy. 
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to a non-negative value, the row marked S, is temporarily chosen as the 
objective row to be maximized. This row is maximized using the normal 
algorithm. However, only the constraints that are currently feasible 
are considered for pivoting. The set of currently feasible constraints 
is referred to as the subproblem. The temporary objective function Sy is 
maximized over this subproblem until Sy reaches a non-negative value. 
At that point Sy. is feasible and its associated row is included in the 
subproblem. Thus once a variable becomes feasible it is not allowed 

to go infeasible at a later iteration. Each infeasible variable is 
treated individually in this manner until all infeasibilities have 

been removed and the subproblem contains the entire constraint set. 

At this point attention is returned to the true objective function 


and the normal algorithm is applied to optimality. 
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CHAPTER III 


COMPUTER IMPLEMENTATION OF THE GONZALES-YOUNG ALGORITHM 


One of the primary objectives of the research reported here 
was to develop a packaged computer program for the solution of AILP 
problems by the G-Y algorithm. This program has been written in 
the procedure-oriented ICETRAN (ICES Fortran) language[1]. However, 
the CDL(Command Definition Language) of ICES[2] has been employed 
to facilitate user interaction with the program. The whole package 
has been incorporated as a component (called OPAILP) of the OPTECH 
subsystem of ICES[3]. OPTECH is the mathematical programming sub- 
System containing several other techniques such as linear programming 
and network flow analysis. 

OPAILP allows the user to specify his problem in much the same 
way as he would normally formulate it. A small set of procedure- 
initiating Commands and their accompanying data are all that are 
required to define the problem to the computer and to initiate the 
solution process. In the absence of input error conditions the 
user can expect to receive his answers in a convenient form. With 
the exception of the time and number of iterations required for 
the solution, no data related to the algorithm is output. 

This chapter is intended fo serve as a user's guide for the 
OPAILP package. It begins with a brief discussion of the consider- 


ations that influenced the design of OPAILP's data structure and 


epee 





data management philosophy. It then specifies the procedures that 

are required in the proper use of the program and briefly describes 

the computer routines that make up the CDB's (Command Data Blocks)[2] 
for command processing. An example problem that encompasses all of 

the OPAILP capabilities 1s presented. Computer source program 

listings and detailed program documentation are contained in Appendix A 


of this paper. 


3.1 Data Structure Considerations 

The dynamic array capability offered by the ICES system allows 
considerable flexibility in problem data management. The way in 
which problem data are structured into dynamic arrays can affect 
the efficiency of data management to a considerable degree. Usually 
the data structure comes under closest scrutiny when it is possible 
for data requirements to exceed the storage capacity of primary 
computer memory. This is so because the inefficiencies encountered 
in referencing poorly structured core-resident data are small com- 
pared to unnecessaril: repetitive transfers of data between primary 
and secondary storage. In general, when one knows in advance the 
order and frequency of data requirements of the algorithm under 
consideration, the design of an efficient data structure is a dis- 
tinct possibility. These requirements are met by the G-Y algorithm. 

In its basic form the algorithm has data requirements that 
are relatively well defined. It consists of the repetitive appli- 


cation of four sequential steps (see Section 2.5). The data that 
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these steps require are either all members of a particular row 


or all members of various columns. 


A. 


D. 


Pricing- requires data elements from the first 
row of the tableau. When a negative coefficient 
is encountered, pricing is halted to see if 
‘breakthrough’ can be achieved. 

Finding Most Binding Constraint- requires the 
zeroth column of the tableau and the pivot 
column that is currently under consideration. 
Entries in column zero and the pivot column are 
referenced sequentially from the top of the 
tableau to the bottom. 

Pivot Acceptance/Rejection- no additional data 
is required. If the pivot column currently 
under consideration generates a zero cut and there 
are more negative columns to be examined, the 
current proposal is rejected. When all negative 
columns have been considered and only zero cuts 
are available, one of these is accepted for 
pivoting. Breakthroughs are always accepted. 
Tableau Updating- requires the referencing of 
all columns of the tableau. The pivot column 

is used to update all other columns by the 


formulae given in Section 2.2. 


When data must be referenced through a dynamic array pointer 
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structure, the number of levels in this structure should be kept to 
a minimum. Thus there is always a question of how the data elements 
should be grouped to form larger referencing entities. For example, 
every step of the G-Y algorithm references data elements as members 
of a particular column or members of a particular row. The 
requirement for one member of a column (row) usually implies an im- 
pending requirement for all other members of the column (row). 
Therefore, grouping the data in the form of row entities or column 
entities seems an obvious possibility. However, there is also the 
possibility of transmitting even larger blocks of data between 
primary and secondary storage (e.g. strings of columns or rows). 

All of these arrangements are possible and the actual steps that 
the algorithm performs provides the basis for meaningful comparisons 
among the possible alternatives. 

The data structure that has been chosen for OPAILP considers 
Single columns as entities in transferring data between primary and 
secondary storage. The idea of row entities was rejected out of hand 
primarily because of Step B presented above. Finding the most binding 
constraint would require the referencing (and transfer) of every 
tableau element each time a particular column was considered for 
pivoting, even though this step requires elements from only two columns 
for this purpose. Since several zero cuts are often rejected before 
a pivot step is actually performed, row-wise referencing would require 
several transfers of every row at each pivot step. It was obvious 


that row-wise referencing was too inefficient to receive serious 
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consideration. This left a choice between single columns and strings 
of columns. This choice was somewhat more difficult, because there 
was a delicate tradeoff involved. Naturally, the larger the entity, 
the fewer will be the number of secondary storage references in any 
given step of the algorithm. However, single column referencing 
provides some types of flexibility that stringing columns does not 
allow. To bring this point to the fore, some more detailed consider- 
ations of the algorithm must be presented. 

There are two aspects of the algorithm that require a consider- 
able degree of flexibility from the data structure. The first of 
these concerns convergence considerations. Both Young's convergence 
proof [4] and the convergence heuristics of Gonzales [5] require the 
logical switching of columns of the tableau. The details of these 
requirements are too elaborate to present here. However, in general 
terms, convergence considerations require that, at some particular 
stage of the computation, column j and column k must logically change 
places in the tableau so that they can be considered for pivot 
selection in reverse order. The second aspect of the algorithm that 
affects the data structure concerns the use of the algorithm as 
a component in the solution of a MILP problem using a method de- 
veloped by Benders [6]. In this context the G-Y algorithm would 
have to add rows to its tableau periodically. The author feels 
that this will eventually be a very important function of the G-Y 
algorithm since all of the requirements of Benders' algorithm are 


currently available in OPTECH except for a required interface 
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between the LP and AILP routines. 

All dynamic array structures are permitted to grow. However, 
if the columns were to be strung together as transfer entities, the 
columns within any particular string would not be able to grow 
conveniently. One would have to either attach unused storage lo- 
cations at the end of every column in anticipation of the rows that 
would be added (see Figure 3.1.1) or provide some way to continue 
each column outside the string (Figure 3.1.2). Either of these 
approaches would be cumbersome; the first because there is no way 
to predict the number of rows that would be added, and the second 
because data management in this structure would be expensive. The 
logical switching of columns in either of these schemes would have 
to be done by physically moving them or by an elaborate bookeeping 
system. On the other hand, single column entities handle both of 
these problems with relative ease. Figure 3.1.3 shows that unlimited 
column growth is permitted, and Figure 3.1.4 demonstrates the ease 
with which two columns can be logically exchanged via the ICETRAN 
SWITCH command [7]. These considerations led to the choice of single 
columns as data transfer entities. However, this choice led to a 
rather detailed data management strategy designed to minimize the 
cost of storing and retrieving data elements in these relatively 


small blocks. 


3.2 Data Management Strategy 


Generally speaking, a good data management strategy will keep 
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the most urgently needed data in primary memory so that accesses 

to secondary storage are kept to a minimum. It will also take into 
account the potential effects of primary memory reorganizations on 

the efficient use of the available data pool [8]. It must predict 

the future need for data of various types and take steps to insure 

that data handling will result in an optimal use of primary memory 

reflecting these needs. 

The ICETRAN language allows the programmer to assume some of 
the data management responsibilities usually performed by the system. 
The programmer is allowed to specify the relative importance of each 
data entity. Each entity may be given a status (released or un- 
released) and a priority (high or low). The data pool may contain 
all types of data with all combinations of status and priority. 
However, when new data are transmitted to primary memory with no 
Space available in the data pool, the monitor system is forced to 
reorganize the pool to create the needed space. Four distinct and 
increasingly drastic options may be exercised by the system for this 


purpose. These four ‘levels of reorganization’ consist of: 


Reorganization Level Action 
0 reposition arrays in pool 
1 purge low-released arrays 
2 purge high-released arrays 
3 purge low-unreleased arrays 


Reorganization at any level purges all arrays effected by that level. 
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The data management strategy that is employed by OPAILP utilizes 
its knowledge of memory reorganization policies in attempting to keep 
urgently needed data in primary memory. The strategy considers that 
there are only two basic operations that the G-Y algorithm perfoms. 
These are pivot selection (Steps A, B, and C) and updating (Step D). 
The updating operations require each of the columns (one at a time). 
This implies that the data can be ‘prepared’ during the updating pro- 
cess for efficient referencing in the pivot selection operation. The 
pivot selection operations consider only negative columns. Other 
columns need never be referenced in this operation. To take advantage 
of this fact OPAILP keeps a copy of the objective function in a one- 
dimensional array outside of the regular tableau structure. This 
array stays at an high-unreleased level! so it can never be removed 
from the data pool]. The algorithm searches this one-dimensional array ~ 
in order to propose pivot .2!umns. Thus a column need never be 
referenced simply to see if it is negative and as such eligible for 
pivot selection. Releasing policies within the updating operations 
insure that there is always a high density of data in the pool. 
Notice that if all columns were released at the same level in the up- 
dating process a memory reorganization would leave the data pool 
completely empty. This would result in grossly inefficient use of 
space. To combat this effect the algorithm keeps track of the number 
of columns that will fit into the pool. It will not release data 


until that limit is approached. However, when the data pool is judged 
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to be relatively full, all other columns that are encountered in 
updating are set to a low-released level. If and when a reorganiza- 
tion is required, it will not occur above Level 1 and space will be 


freed though the pool will remain relatively full. 


3.3 OPTECH/OPAILP Commands and Their Uses 


The OPTECH executive routine provides the set of command data 
blocks used by the ICES executive program in processing OPTECH com- 
mands. This executive routine is written in the CDL language of 
ICES. The OPTECH/OPAILP component encompasses only six of these 
commands. A user-constructed set of commands, and command-dependent 
input data, is used to initiate the OPTECH/OPAILP routines. The set 
of commands listed on the following page, and their associated data, 
constitute the full capability of OPAILP. The Required commands are 
the minimal set that can be used to solve an AILP problem. They 
are presented in the order that they must occur in the user's data 
deck. Any other ordering of these commands will cause a termi- 
nation of the job. The Optional commands provide additional capa- 
bilities. They also have a required ordering, but this will be 
discussed in the next section where sequencing and data requirements 


are specified. 


Oe 


COMMAND 


FUNCTION 


REQUIRED COMMANDS 


OPTECH 
AILP *ROW ny ‘COL No" ‘LIMIT he 


RIGHT SIDE 
{rhs data deck] 
TABLEAU 


[Matrix data deck] 
ITERATE 


FINISH 


initiates OPTECH subsystem 


initiates the G-Y procedure 

The terms in '. ‘ are optional. 
They specify the number of rows, 
the number of columns, and the 
user's limit on the number of 
iterations that will be performed. 


initiates the input of rhs data 


initiates the input of the matrix 
coefficients 


initiates solution process 
Standard output occurs when the 
process is completed. 


system command marking the end 
of the job 


OPTIONAL COMMANDS 


INSERT ROW ny 


[New row data deck ] 


WRITE (Any OPAILP command) 


initiates the addition of a row 
The row is added to an optimal 
tableau and a reoptimization is 
performed. 


echo=-prints the data accompany- 
ing any OPAILP command 





3.4 Command Documentation 
This section explains the function of each of the OPAILP 
commands. The underlined portion of each command name constitutes 
a short form version of the command that is acceptable. All data 
are read in under fixed format conditions by the CDB's associated 
with each command. This restriction has been imposed because the 
volume of data that is usually required would lead to prohibitive 
input times if the data were to be processed under the free format 
that is available with the CDL. Particular consideration should be 
given to these format restrictions since all data is input in either 
integer or alphanumeric form. In every case the data should be 
right justified in the data field that is provided. 
Command Documentation OPTECH: 
1. External command structure: 
OPTECH 
2. Description: The system command OPTECH must be the 
first command in any OPAILP run; it requests the ICES 
system to transfer control to the OPTECH subsystem 
(of which OPAILP is a part). 
Command Documentation AILP: 
1. External command structure: 


AILP ROWn, COL n, LIMIT ng 


1 


2. Description: The AILP command immediately follows the 
OPTECH system-level command; it initiates the G-Y AILP 
procedures. 
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The optional modifiers are used as follows: ROW ny Inputs the 
user's estimate of the number of rows, Ny° COL No inputs the number 
of columns, No LIMIT N3 inputs a limit on the maximum number of 
iterations that the user will permit. If the problem is not solved 
within this limit, the best available solution will be output. 

Though the problem size modifiers are optional, their omission 
May increase the solution time slightly. The numbers n and Ny are 
meant to measure the size of the problem as formulated (and not to 
estimate the size of the tableau). They are used to initialize the 
dynamic array structures. In the absence of user specification, 
Standard values of 5 are assumed. Use of the iteration limit (n3) 
Should be considered. No matter how small the problem, there is no 
way to estimate the number of iterations that will be required for 
an optimal solution. A standard value of 100 is assumed when the 
user does not specifically set Ng 
Command Documentation RIGHT SIDE: 

1. External command structure: 

RIGHT SIDE 

2. Description: The RIGHT SIDE command initiates the 

input of right hand side data. These are the b coef- 
ficients in the problem formulation of Section 2.1. 
These data are placed in the zeroth column of the 
tableau and the row names that accompany each of the 
coefficients are used to guide the input of matrix 
data under the TABLEAU command. One rhs coefficient 


(whether zero or not) must be input for each constraint 
row. The objective is not input under this command. 
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3. Format for rhs data deck: Three pieces of data 


Command 


1 


are required for each constraint. These three 
data are: constraint type (L for <, G for >, 
or a blank field to indicate an equation); row 
name (a unique name of up to four characters 
for each row); and the coefficient itself (an 
integer number). The L or G of inequality type 
must appear in card column 8. The row name must 
appear within the field between card columns 9 
and 12. The integer coefficient must be right 
justified to card column 20. The occurrence of 
any character in card columns 1 through 4 will 
mark the end of rhs data (e.g., END). Refer to 
the example problem of Section 3.5. 


Documentation TABLEAU: 
External command structure: 
TABLEAU 


Description: The TABLEAU command initiates the 
input of matrix coefficients data. These are the 


as of the problem formulation in Section 2.1. 


There is a variable name and a row name associated 
with every matrix coefficient. The row names are 
matched against row names input.with the RIGHT SIDE 
command in order to position the coefficients in 
the proper location in.the tableau. The variable 
name defines the names of the problem variables. 
Tableau input is column dependent and all coeffi- 
cients that are associated with a particular 
variable must be contiguous-.in the tableau data 
deck (though row names within.any column group may 
occur in any order). 





Format for Matrix data deck: . The variable name can 
be up to four characters and must.appear within the 
field between card columns 5 and 8. JHowever, the 
variable name must be positioned in that field in 
exactly the same manner every.time. The row name 
associated with the tableau.coefficient must appear 
in the field between card columns 9.and 12 and must 
appear in exactly the same manner. as it did in the 
RIGHT SIDE data deck.. The row. name associated with 
the objective function is always OBJ and it must be 
right justified to column 12. The coefficient itself 


02 





Command 


a) 
° 


Command 


1 


a 


must be an integer number right justified to 

card column 20. Matrix coefficients that have 
zero values need not, and should not, be in- 
cluded in the TABLEAU data deck. Slack variables 
are automatically generated for each row. The 
user should not instal] slacks before inputting 
the constraints since this will lead to a Phase I! 
being required (equations). The occurrence of 
any character in columns 1] through 4 will termi- 
nate the input process (e.g., END). 


Documentation ITERATE : 
External command structure: 
ITERATE 


Description: The ITERATE command simply initiates 
the solution process. Provided the RIGHT SIDE and 
[TABLEAU commands have occurred in proper sequence, 
the ITERATE command will attempt to maximize the 
objective function. Output will immediately follow 
this command whether the optimal solution is 
achieved or the user's iteration limit is exceeded. 


Documentation INSERT ROW: 
External command structure: 
INSERT ROW 'TYPE®* 


Description: The INSERT ROW command may occur at any 
time after an optimal solution is attained. (Several 
rows may be added in succession.) Only one row may 

be added with each INSERT ROW command. This command 
takes a constraint expressed in terms of the original 
problem variables and updates it to the optimal tableau 
representation. If the insertion of a row causes an 
infeasibility to appear, the tableau will be re- 
optimized and the new solution_.will be output. The 
term 'TYPE' refers to the type.of constraint that is 
being added (e.g., L, G, etc.). The TYPE data is 
enclosed in single quotes. and. must be four characters 
in length (i.e., just as.the. constraint type was speci- 
fied in the RIGHT SIDE input). 
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3. Format for INSERT ROW data deck: The coefficients 
of a new row may be read in any order. There 
1S a row name and column name associated with 
every coefficient. The data format for INSERT 
ROW data is exactly like the data format for the 
TABLEAU command. However, the rhs value for a 
row that is betng inserted must be given the 
column name OBJ right justified to card column 8. 
Every row coefficient (whether zero or non-zero) 
must be input explicitly and must be right justi- 
fied to card column 20. Any character appearing 
tn card columns | through 4 will terminate the 
INSERT ROW input data (e.g., END). 


Command Documentation WRITE: 
1. Externa! structure: 
WRITE (RIGHT SIDE, OR TABLEAU, OR INSERT ROW) 
2. Description: The WRITE command may precede any 
of the commands listed in the parentheses in 


order to cause the echo-printing of their data 
decks (e.g., WRITE TABLEAU). 


3.5 An Example Problem 


In order to provide a better insight into the procedures required 
in using OPAILP, an illustrative problem will be presented. The 
function of each command that is used in this problem will be ex- 


plained. The problem under consideration is the following 
max Ax, + 3X. + 2X4 


SUG3: && eK) + 3X. - X3 ¢ 9 
- 3x, € 2x5 + eX = 12 


4x, = 3x, * 


Pan = 


le 





Notice in particular that the second constraint is an equation so that 
Phase I is required. As a point of interest it should be noted that 
the pivot steps employed in gaining a feasible solution to this equa- 
tion result in an infeasible solution to the third constraint requiring 
the other part of Phase I described in Section 2.6 to be employed. 
Figure 3.5.1 is a facsimile of the computer printout that will result 
from proper specification of the problem mentioned above, The function 
of each of these pieces of data will be explained. The user commands 
are assumed to have begun in card column 1, though there are no format 
restrictions. 

The first command, OPTECH, causes a transfer of contro! to the 
OPTECH subsystem. The AILP command indicates a 3 x 3 constraint 
matrix and a user imposed limit of 20 iterations. WRITE RIGHT SIDE 
causes the echo-printing of the rhs data deck. In this data deck 
notice that the second constraint is specified as an equation by the 
absence of an L or G in column 8. The constraint row names are 1, 2 
and 3. WRITE TABLEAU causes the echo-printing of the matrix data deck. 
The row name of the objective function is always OBJ while the other 
row names match the row names of the rhs data deck. The ITERATE com- 
mand causes the solution procedure to be initiated. The optima! solu- 
tion follows. Notice that slacks have been insta led and numbered in 


the order the constraints were input. 
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The INSERT ROW conmand initiates the addition of a new 


inequality: 


xy + Xo + X32 < 14 


Note that OBJ is used in the data deck to indicate the column name 
of the rhs value of this inequality. Since the previous solution 
does not satisfy this constraint, a reoptimization results and the 
new optimal solution follows. The fourth problem constraint has 
been added to the bottom of the tableau and its slack S004 has been 


generated. 





OPTECH 


AILP ROW 3 COL 3 LIMIT 20 


WRITE RIGHT SIDE 


RIGHT SIDE 
Lo 
2 
ie 
WRITE TABLEAU 
TABLEAU 
Xl OBJ 
Xl 
Xl 2 
Pe a3 
X2 OBJ 
 — 
KZ. 2 
XZ 
X3 OBJ 
koe 
od 
x3. 3 
ITERATE 


FEASIBLE ITERATION 


12 
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TOTAL ITERATIONS 
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OBJ 4 
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3.6 Output Messages 

The following is a list of all output messages that OPAILP may 
print and the situations under which they occur. 

1. FEASIBLE AT ITERATION n 


Printed when the algorithm has attained an 
initial feasible solution 


2. OPTIMALITY HAS BEEN PROVED 


Printed immediately when an optimal tableau is 
achieved 


3. TOTAL ITERATIONS ny 


Printed before output of problem variables 
4. STATIONARY ITERATIONS no 

Printed before output of problem variables 
5. USER ITERATION LIMIT EXCEEDED 


Printed when optimality has not been achieved 
within the user-supplied iteration limit (or 
100 iterations if the user does not specify a 
limit) 


6. COMMAND OUT OF SEQUENCE 


Printed when a command is encountered before 
jt can be processed (Termination always follows.) 


7. NO MATCH FOR — name — IN ABOVE CARD 


Printed whenever a column or row name in the TABLEAU 
or INSERT ROW data decks cannot be matched with a 
name that has already been input. (This command 

1s always preceded by an echo-print of the card in 
which the error was found. ) 


8. TOO FEW (MANY) COEFF. IN INSERT ROW COMMAND 


Printed when INSERT ROW data deck does not supply 
a row that is compatible with the problem 
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12. 


JEP 


14, 


IN ABOVE CARD, ROW (COL) NAME IS REQ. 


Printed when rhs or tableau data card does not 
Supply a row (column) name when needed 


- IN ABOVE CARD — — IS ILLEGAL 


Printed when a rhs data card does not supply 
an acceptable inequality type(Card precedes 
the message. ) 


» NO INTEGER SOLUTION EXISTS TO EQUATION name 


Printed when Phase I fails to find an integer 
solution to the equation named on the message 


THIS PROBLEM IS UNBOUNDED 


Printed when it is determined that there 1s no 
upper bound on the value of the objective 
function (Job is terminated with no output.) 


THIS PROBLEM HAS NO FEASIBLE SOLUTION 


Printed when Phase I is unable to raise a 
negative slack to a non-negative level (Job is 
terminated. ) 


JOB FLUSHED 


Printed with most error conditions to indicate 
that the entire job has.been scrubbed and EXIT 
has been called. 


Any input error that violates a format restriction will result 


in a FORTRAN error code and control will be taken away from OPAILP. 


The entire job will be terminated. 


3.7 Some Tips on the Use of OPAILP 
The G-Y algorithm has particular difficulty with a certain type 


of problem. The inclusion of a constraint row that contains coefficients 
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that are rather large and are relatively prime to each other will 
almost invariably lead to a large number of iterations. Jhe reasons 
behind this characteristic would require a rather detailed explana- 
tion that probably would not interest the user of this subsystem. 
However, the user should bear this point in mind when formulating 
his problem. The integer coefficients should be made as smal] as 
possible, commensurate with the accuracy of modeling that is re- 
quired. 

The inclusion of equations in the constraint set often produces 
a problem that is somewhat more difficult to solve than problems 
with only inequalities. Thus the user should never add slack vari- 
ables to inequality constraints before they are input to the program. 
This will initiate a rather expensive (and unnecessary) Phase I 
attempt to get a starting basic feasible solution. 

If the AILP command includes an accurate estimate of the number 
of rows and columns of the problem, some computation time will be 
saved. Estimates that are too low will cause the redefinition of 
dynamic array sizes and this takes additional time. 

The WRITE command allows the user to receive output that shows 
both his problem definition and his results. However, this feature 
Should not be used if its purpose is to check for keypunching errors 


in the alphanumeric fields. When keypunching mistakes lead to errors 
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in these fields the program always echo-prints the particular card 
that contains the error along with a message identifying the mis- 
take. Input-output operations are always expensive and should be 


avoided where possible. 
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CHAPTER IV 


CONVERGENCE DIFFICULTIES IN THE GONZALES-YOUNG ALGORITHM 


This chapter presents a geometric interpretation of the G-Y al- 
gorithm and isolates various phenomena that cause its computations to 
be unnecessarily prolonged. These observations are subsequently used 
in Chapter V to present a revised concept for a primal AILP algorithm 
and to develop an algorithm which specifically avoids certain in- 


efficiencies occurring in the G-Y algorithm. 


4.1 Convex and Restricted Polyhedra 


It is a well known fact that associated with any set of linear 
constraints there is a unique ‘convex polyhedron’, the vertices of 
which correspond to the basic feasible solutions to the associated 
linear program. The 'faces' of this polyhedron are defined by the 
geometric representation of the problem constraints as hyperplanes 
in n-dimensional space. The intersection of n-1] hyperplanes in n- 
dimensional space defines a straight line segment. A line segment 
that contains two vertices of the convex polyhedron is called an 
‘edge' of that polyhedron. Dantzig [1] has offered a geometric 
interpretation of his Simplex method for linear programs in terms 
of n-dimensional non-basic (parametric) space. This interpretation 
views a Simplex pivot step as translating the non-basic coordinate 


system from one vertex of the convex polyhedron to an adjacent 
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vertex by moving it along the edge of the convex polyhedron 
connecting the two vertices in question. This translation is 
accompanied by a rotation of the coordinate axes so that the 
geometric representation of the convex polyhedron is distorted 
as a result. This geometric interpretation was displayed in 
Section 2.3. 

When the same set of linear constraints is used in the 
formulation of an AILP problem the vertices (basic feasible 
solutions) of the convex polyhedron may no longer be of interest 
because they may not be integer points. However, all of the 
integer points that satisfy the constraint set must be contained 
in the convex polyhedron. Furthermore, associated with any con- 
vex polyhedron containing a set of integer points there can be 
only one ‘restricted polyhedron’ having the property that every 
vertex is an integer point from the set and the entire set is 
contained within. The concept of a restricted polyhedron is a 
useful concept in a geometric interpretation of the G-Y al- 
gorithm. Clearly if all of the constraints corresponding to 
the hyperplane faces of the restricted polyhedron could be 
added to the original constraint set a normal application of 
the Simplex method would result in only integer solutions to 
the problem. However, the G-Y algorithm recognizes that there 
are other, less complicated, ways of attaining integer solutions. 
Rather than attempting to generate actual faces of the restricted 


polyhedron, the algorithm simply attempts to uncover edges. There 
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is a whole family of hyperplanes in n-dimensional space (n > 2) 
that contain a given edge of the restricted polyhedron and an in- 
finite number of hyperplanes from this family border closed hal f- 
Spaces that contain all of the feasible integer points. Pivoting 
on any n-1] of these hyperplanes in succession will ‘uncover’ a 
given edge of the restricted polyhedron. This is not to say that 
the algorithm has a priori information about the restricted poly- 
hedron. It only means that the scope of the task that the cuts 
must perform is not as broad as might be anticipated. The G-Y 
algorithm, like all other cutting plane algorithms, is merely a 
modified Simplex method. The cuts that are generated function 
primarily to surpress non-integer solutions that the Simplex 

method would normally produce. Actually the term ‘cut’ appears 

to be something of a misnomer because these additional constraints 
do not always cut away portions of the convex polyhedron (this 
topic will be pursued in the succeeding paragraphs). However, it 
is important to realize that these cuts are the only thing that 
distinguishes the G-Y algorithm from the Simplex method. Therefore, 
any difference in efficiency between the G-Y and the Simplex method 
1S purely a function of the effectiveness of the Gomory cuts that 


are employed in the G-Y algorithm. 


4.2 Convergence Difficulties - An Interpretation 


As yet no method has been devised by which properties of the 


restricted polyhedron can be extracted from computations involving 
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the original constraint set that defines the convex polyhedron. 
However, it 1s possible that in some cases the restricted poly- 
hedron of a given constraint set might be very similar to its 
convex polyhedron. In fact there might be cases in which the 
two polyhedra coincide. One would surmise that these cases would 
yield quite readily to the Simplex-like G-Y algorithm. It seems 
reasonable to assume that the number of iterations required to 
solve the AILP problem would compare favorably with the number 
of iterations required by a Simplex solution of the associated 
LP problem. Investigations based on this premise have shown 
this assumption to be false. Several small problems were 
formulated in such a way that the convex and restricted poly- 
hedra were identical. On these problems the G-Y algorithm 
usually required many more iterations than the Simplex method 
even though both methods yielded exactly the same integer 
solutions. These investigations led the author to recognize 
several important aspect of the G-Y algorithm, each of which 
contributes to unnecessary inefficiency. The first of these 
was a recognition that the G-Y algorithm often generates com- 
pletely superfluous cuts. This topic is taken up in the remainder 
of the current section. Discussion of the rest of these points 
will be postponed until Chapter V where they are used to develop 
a modified algorithm. 

One of the most striking things that was observed in the 


G-Y solution attempts on the problems mentioned above was the 
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fact that a very high percentage of the cuts that were generated 

were completely redundant. In other words these cuts excluded no 
portion of the convex polyhedron. Perhaps more interesting was the 
fact that these redundant cuts tended to follow one another in 
succession. In other words the occurrence of one redundant cut was 
most often followed by a whole string of cuts of the same type. The 
hyperplanes associated with the redundant cuts in such a string al] 
intersected the (modified) convex polyhedron at the same integer 
vertex. All of these cuts (except perhaps the first) were zero cuts. 
Investigations have shown that these redundant cuts can be anticipated 
and, with major revisions to the algorithm, they can be avoided. 
Because the geometric interpretation will be so important in under- 
Standing the reasoning behind the modifications that will be proposed 
later, it will be helpful to first examine the effect that this tyne of 
cut can have on the rate of convergence of the G-Y algorithm before 
determining how they can be anticipated. 

Recall that the intersection of n hyperplanes in n-dimensional 
Space defines a point in that space. If only n hyperplanes associated 
with the problem constraints intersect at each vertex of the convex 
polyhedron no special difficulty arises in solving the LP problem. 

In this case there ére (=) = n edges of the convex polyhedron 
eminating from each vertex. Pivoting in any one of the n columns of 
the tableau will translate the non-basic (parametric) coordinate 
System down one of these edges to an adjacent vertex. Pivot selection 


rules are designed in such a way as to improve the objective function 
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at every iteration. Since there is only a finite number of vertices 
there can be only a finite number of iterations before optimality is 
attained. On the other hand when more than n hyperplanes, say ntq, 
intersect at a given vertex there are (Ty) line segments eminating 
from that vertex. However, not all of these line segments can be 
edges of the convex polyhedron. In fact only ntq-1 of these line 
segments can be edges of that 2): hecron. The other line segments 
lead from the vertex in question out into the infeasible region and 
an attempt to translate the non-basic (parametric) coordinate system 
along one of these line segments will be unsuccessful. Proper 
application of the Simplex pivot rules will prevent such movement 
and a degenerate solution will result (a change in basis with no 
change in basic feasible solution). Geometrically speaking Dantzig's 
lexico-graphic pivot selection rules monitor the attempted movements 
down the various line segments. They require that the line segments 
be explored systematically so that no particular movement is ever 
attempted twice. Since the number of these line segments is finite, 
eventually one that is an edge of the convex polyhedron must be ex- 
plored. When this occurs translation of the coordinate system is 
achieved and again the algorithm is guaranteed to be finite. 

Applying the same geometric interpretation to the G-Y algorithm 
immediately points to the convergence problem. Consider the possibility 
of redundant cuts of the nature previously described. All of these 


cuts intersect at the same integer vertex. At most one of them was 


needed by the algorithm to define that vertex (only the first of them 
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was generated in pivoting to the vertex in question). Furthermore, 
none of these redundant cuts can help to define the next edge of the 
restricted polyhedron so pivoting on these cuts only serves to per- 
petuate the degeneracy already existing at that integer vertex. The 
Simplex lexico-graphic pivot s:lection rules no longer suffice since 
each new redundant cut adds new line segments at the vertex, none 
of which are edges of the restricted polyhedron. 

Though Young [ 2] did not address himself to the possibility 
of redundant cuts in the algoritm he did show that the algorithm 
could be made to converge in a finite number of iterations. However, 
it is the rate of convergence of the algorithm and not its ultimate 
convergence that is of primary interest here. With this in mind it 
will first be necessary to establish how these redundant cuts arise 
in the algorithm before methods can be proposed to avoid them. The 
following section takes up the first of these issues; the second is 


the topic of Chapter V. 


4.3 Occurrence of Redundant Cuts in the Gonzales-Young Algorithm 


Recall that each of the problem variables may be viewed as 
Slacks associated with hyperplanes in the n-dimensional non-basic 
(parametric) space. Each constraint row in the tableau represents 
Such an hyperplane. A constraint (and hence an hyperplane) will be 
referred to by the name of its slack variable. Points that satisfy 
a given constraint yield non-negative values in its slack variable. 


A constraint S. 1s redundant to another constraint Si. when al] 
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non-negative points that satisfy Sy. also satisfy Si > but non- 
negative points that satisfy S ; do not necessarily satisfy Sp 
The phrase ‘do not necessarily satisfy...’ is used because it wil] 
allow S. to be redundant to Si. if it is a scalar multiple of Spe 
In avoiding redundant cuts it will also be desirable to avoid cuts 
that duplicate constraints that are already present in the tableau. 
Returning to the notational conventions of Chapter II, the 
following statement prescribes the circumstances under which redundant 
cuts are generated by the G-Y algorithm: 
“Any cut that 1s generated in column k from row 7 is redundant 
if ai 9/45 k = integer" 
Proof: 


Consider the constraint: 


si, as 9 + ay j(-t,) (4332) 
j=! 
and suppose that a Gomory cut 7s generated from S; in 
column k (i.e. pas in (2.3.2)). That cut would be: 
n 
r= Lay g/ay lt) Cay j/ay yl (-t,) (4.3.2) 
j=l 


By definition of la; 5/4; J , the as 120, lzcee6n 


in (4.3.1) can be rewritten as: 


455 7 84k bE f/4G kd * FG, 


where: 0 < fig & ask 


Jaye 





Thus (4.3.1) can be rewritten as: 


n 


* tumeeieisk {C35 0/244 +) Lay 5/tiall-ty)} 
j=] 
n 
+ Fi 0 + > fF; j(-t,) (453.3) 


ie. 


The following observations can be made about (4.3.3): 
(1) the term inf fis equal to r 
(ii) a, 9/4; x = integer, implies fio = Q 
(iii) a; , > 0 by assumption 


therefore (4.3.3) can be written as: 


n 


S$, = a; ,(r) + > Ti 4(-t.) (4.3.4) 
j=] ’ 


Consider the t. > 0 j=l,...,n that satisfy the 

cut r (i.e. r > 0). For such points S; is not necessarily 
non-negative. In fact if S. is not an integer multiple of 
r (i.e. f; . = 0 Y¥j) then there are points lying on the 


T5J 


hyperplane r (r=0) that do not satisfy s, (i.e. s; < Oe 


Now rewrite (4.3.4) in terms of r: 


n 
= (1/a, Tey 3. : f; jes) (4.3.5) 
j= 


And consider E > 0 j=l,...,n that satisfy the constraint Ss. 
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Such points yield S; > QO and since all other terms on the right 
side of (4.3.5) are non-negative, they all satisfy the constraint r 
as well. Therefore r is redundant to Si. 

Recalling the previous discussion concerning convergence 
difficulties, the effect of generating and pivoting on the cut r 
can now be explored. Suppose that the tableau is currently all- 
integer. Let the superscript i on entries Gs imply that ie 
is an entry in the ie tableau and superscript i+] imply tableau 
i+], etc. By assumption that as o/s = integer , it is obvious 
that as /4; = Lay 9/8; yd. Futher observe that pivoting on 


column k in tableau 7 will result in the following situation in 


tableau i+]: 


r= anneto = £2) 0/84 q) ~ C25 4/9 IL; 9/8} gl = 0 
(refer to updating formulae in Section 2.2) 

Thus the cut r is not only redundant to the constraint Sis but 
both S; and r intersect the same vertex of the convex polyhedron 
(the vertex at the origin of the non-basic coordinate system of 
tableau i+1). Therefore, pivoting on the cut r in tableau i has 
added degenerate solutions to the problem in the manner discussed 
in the previous section. 


Having observed the circumstances under which redundant cuts 


are generated, attention can be turned to the somewhat more interesting 
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topic of redundant cuts occurring in strings. It seems reasonable 
to assume that pivoting on one redundant cut somehow leads to the 
generation of additional redundant cuts in succeeding tableaus. 
This assumption is correct and the circumstances that lead to this 
result can be stated as follows: 

"If a cut is generated in column k from row i having the 


> 


property a. 0/4; a integer and there are entries a. 


in the row such that a fa: k # integer; each such 


15] 
column j will generate a redundant cut from row 1 in 
the succeeding tableau." 
Proof: 
Consider the tableau updating formulae in Section 2.2. In 


row i of the updated tableau (it1): 


re . . . 
ae = ay 9 - a} lay o/ay gd = 0 (4.3.6) 
ye it] _ oi jr ae 4.3.7 
and: a 45 a; lay 3/4; 4! (4.3.7) 
, 1 im _ A j j 
i 45/8 yk = bay 3/ay gl tf; 
; i | 
-° arsine J 
And: eas = 0 (integer), and by the previous proof 


column j generates a redundant cut from row i in tableau i+]. 
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Thus it has been seen that redundant cuts rot only create 
degenerate solutions, but also can lead directly to the generation 
of still more redundant cuts and degeneracy is perpetuated through 
the algorithm. With this thought in mind one can easily appreciate 
the magnitude of the effort represented in Young's proof of finite 
convergence. However, it is still reasonable to ask why the generation 
of redundant cuts is tolerated in the algorithm and to ask whether 
there might be some way to avoid them. 

Suppose that at some stage of the algorithm a constraint Ss; is 
encountered such that S. would generate a redundant cut r. It has 
already been seen that both r and Ss; intersect a common vertex and 
r is redundant to Ss It is logical to ask why the algorithm could 
not simply neglect to generate r and pivot on S, instead. The reason 
that this alternative is not taken is that ask (the pivot element) 
is not necessarily equal to +] and though the updated tableau would 
have an integer solution, the other tableau entries would not be 
integers. Therefore, succeeding pivot steps could not be guaranteed 
to produce integer solutions. At that point the algorithm (in its 
present form) would break down. 

The possibility of avoiding redundant cuts in the manner described 
above led the author to search for methods for restoring the integrality 
of the tableau once it has been lost. Such a method was discovered, 
but it was found to generate redundant cuts much in the same way that 
was described above, so it was abandoned. However, the concept of 


this method was derived directly from the author's realization that 
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the G-Y algorithm is rather strongly related to Gomory's Algorithm I 
(a dual feasible method) though on the surface the two algorithms 
seem completely dissimilar. This observation has led the author 

to propose a modified primal AILP algorithm that draws from both 

of the algorithms mentioned above. This observed relationship 
between the two algorithms will be presented first in Chapter V 


before the modified algorithm is described. 
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CHAPTER V 


A MODIFIED PRIMAL AILP ALGORITHM 





This chapter draws on the previous geometric interpretation 
of the G-Y algorithm in order to motivate the reasoning behind the 
modifications that are proposed. It was seen in Chapter IV that the 
G-Y algorithm tolerates the generation of redundant cuts in order 
to preserve an all-integer tableau. The present chapter shows that 
an all-integer tableau is not essential to the concept of a primal 
AILP algorithm. Methods are presented here that not only handle 
the non-integer tableau, but also specifically avoid the generation 
of redundant cuts and their accompanying inefficiencies. Gomory 
cuts remain the vehicle by which integer solutions are produced, so 
essentially no new concepts are introduced in the context of Integer 
Linear Programming. However, since the cuts that are used in the 
modified algorithm are developed in a manner somewhat different from 
the way cuts are developed in the G-Y algorithm, it will be necessary 
first to present these cuts in their original context (Gomory's 
Algorithm I). Then, in order to bring the geometric interpretation 
abreast of the theoretical development, this chapter shows how the 
different cuts can be used to produce a primal algorithm that is 
equivalent to the G-Y algorithm. This chapter goes on to generalize 
the concept of this equivalent algorithm to produce a modified primal 


AILP algorithm. This modified algorithm is compared to various other 
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AILP algorithms (including G-Y) on several problems. Finally an 
attempt is made to explain the different rates of convergence among 
the algorithms and to place the modified algorithm in perspective 


with respect to other cutting plane methods. 


5.1 Review of Gomory's Algorithm I 
Chapter I traced the development of cutting plane methods to 


the present. Recall that Gomory's Algorithm I was the first proven 

algorithm for the AILP problem. This algorithm is a dual feasible 

method in that no integer solution is produced until the optimum 

is found. The integer programming portion of the algorithm begins 

with an optimal tableau for the associated LP problem. The generation 

of each cut is followed by a complete reoptimization of the newly de- 

fined (by virtue of the new constraint) LP problem. Optimality is 

achieved when the reoptimized LP problem has an all-integer solution. 
From the optimal tableau of the LP problem Gomory selected a 

row that represented one of the integer constrained problem variables 


whose value was non-integer. 


n 


HT a gt) ae G(-ty) (5.1.1) 


j=l 
and showed that non-negative integers t, j=l,.c.,n tat 


satisfy (5.1.1) also satisfy: 


read g * Dd 5(-t3) (5.1.2) 





where: dy 5 = a5 7 Cay 5 (2c 


and [a, ; denotes the largest integer _ ayy 


(a, j typically non-integer) 


He also showed that r was constrained to cake on non-negative values 
by virtue of the way it was derived. When (5,..2) was appended to 
the tableau it made the previous optimal solution to the LP in- 
feasible and a reoptimization was required. Geometrically speaking 
the cut (5.1.2) actually cut the previous optimal solution away 

from the convex polyhedron in the same way that the cut (2.3.2) is 
expected to function in the G-Y algorithm. However, it 1s evident 
that (5.1.2) could never be a redundant cut since deo could never 
be zero. The vertex of the convex polyhedron that corresponds to 
the basic feasible solution of the tableau from which r is generated 


is always excluded by the inclusion of r in the constraint set, 


9-2 On the Relationship Between the G-Y Algorithm and Algorithm I 


Since the two types of cuts that have been discussed serve 
roughly the same purpose in the algorithms that generate them it 
seems reasonable to assume that they are more strongly related 
than their derivations tend to indicate. The discussion that 
follows will show that this assumption is correct. 

Recall that discussions in Chapter IV reached the conclusion 


that the G-Y algorithm must tolerate redundant cuts in order to 
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preserve integrality of the tableau. It was seen there that these 
redundant cuts could be anticipated and could be avoided if methods 
were available to restore the tableau's integrality. In searching 
for such a method the author discovered that cuts used in Gomory's 
Algorithm I (5.1.2) could be helpful in this effort. It was also 
observed that these cuts could be used in a primal context to 
produce an algorithm that is equivalent to the G-Y algorithm. This 
equivalent algorithm starts with an all-integer’tableau just as the 
G-Y algorithm. However, it requires two distinct pivot steps to 
achieve the same results that are achieved in the G-Y algorithm in 
one pivot step. The equivalent algorithm first performs a normal 
Simplex pivot step in the first tableau. Then from the updated 
tableau a cut (5.1.2) is generated in a special way. This cut turns 
out to be exactly the same cut as the one that would have been 
generated by the G-Y algorithm in the initial tableau. Therefore, 
pivoting on this cut in the updated tableau produces exactly the 
Same results aS one pivot step of the G-Y algorithm. This procedure 
will be presented in detail] below, but first one important point 
should be considered. Although the equivalent algorithm requires 
two pivot steps for every one required by the G-Y algorithm, its 
procedure does demonstrate an ability to handle a non-integer 
tableau. Thus the equivalent algorithm is presented here only to 
motivate the concepts that will be generalized later to produce a 


modified primal AILP algorithm. 


ve 





Suppose that an initial all integer tableau is given and that 
the G-Y algorithm has decided to generate a cut (2.3.2) in column k 


from the constraint: 


n 


cN 
S$; = aj 9? a a; j(-t,) (5.2.1) 
j=l 
that cut would be 
n 
5 = [ay o/a; 54 + » are Cone (-t,) (5.2.2) 
j=l 


However, suppose that instead of generating and pivoting on (5.2.2) 
a normal Simplex pivot step was performed by selecting ai k as the 
pivot element. This would raise the parametric variable ty to a 
positive level (bring it into the basis) and replace it in the non- 
basic (parametric) variable set with t*. Observe the effect of this 
pivot step on the tableau representation of the variables t, and S.. 


In the initial tableau: 


t, = 0 + (-1)(-t,) (Saco) 
n 

$; = 4; 9 + 2 a. wglots) (5-24) 
j=l 





In the updated tableau: 


n 

th = a; 9/45, * > ag 5/44 glo ty) + May (tp) (5.2.5) 
ms 
jFk 

s. = 0+ (-1) (-tf) (5.2.6) 


Now if a cut (5.1.2) is generated from the updated tableau repre- 


sentation of ty that cut would be: 


n 
! — _ =< = -= -_ 
r' = d; ot) ds 5) t.) + ( Way dK t*) (5.2.7) 
j=l 
tk 
where di = ay 4/4i k - la; j/a; J (5.250) 


We wish to show that the cut (5.2.7) is exactly the same as the cut 
(5.2.2) (i.e., r'=r). Replace the d. j in (5.2.7) with the right 
side of (5.2.8). From (5.2.6) observe that t* = s, and substitute 


the right side of (5.2.4) for tf into (5.2.7) to yield: 


n 
r'=[a; o/a; , J+ ) Lay 5/4 gl (-ty) “Ay g/8y , meee 


j=l 
jFk 
n n 
oo” a, j/a; ,(-t,) +1735 85 oD setp} (5.2.9) 
jel cat 
J#k 
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Cancellation of terms in (5.2.9) yields: 


n 
r's[a; o/a, Jt ) bay, 5/94 gll-ty) Bia ei cea (5.210) 
j=l 
j7k 
and since (a; ./a; ) = La; ,/aj yd = 1, (532510) osimpliitiessto 
r' = r. Furthermore, pivoting in column k on the cut (5.2.7) will 


produce the following tableau representations: 


n 
tg 0/84 (744.0 Via, if tink Sig) ae 
a4 
Aen ae" (5.2.11) 
n \ 
S. = as ,(d i, 9) + as k Y d; jf-t,) - a. ace (Sez, lie} 
a 
J#k 
or by substituting the right side of (5.2.8) for d. ; 
n 
c= fa; 9/45 4d + Q [a, 3, ,J(-t 3) + (-t,) (5 cah3) 
ie 
sk 


~” 
i 
my 


47 9400 7 84 Kb4G 0/84 4! t- 


ost (a; ja; fay 5/a; pd) (-ty) Tal et (5.2.14) 


— 


\f 
But recall from Chapter II that: 


Cas 


45657 84k Fy i/AG al tf, 


ie 





So (5.2.14) can be rewritten as: 


n 
s. =f. ,¢ ) fy j(-ty) + (-t,) (5.2.15) 


By observing the pivot formulae (2.2.2, 2.2.3) one can see that 
(5.2.13) and (5.2.15) are exactly the same results that would have 
been attained by pivoting on the cut proposed by the G-Y algorithm 
in the initial tableau. Thus it is seen that, by using the cuts 
that are generated in Gomory's Algorithm I in a special way, a 
primal algorithm can be produced that is equivalent to the G-Y 
algorithm. This equivalent algorithm works with a non-integer 
tableau in the second step of its procedure. As an aid to under- 
Standing the somewhat cumbersome development presented above, Figures 
5.2.1 and 5.2.2 illustrate the parallel between the two techniques. 
The two tableaus of Figure 5.2.1 reproduce the G-Y pivot step pre- 
sented in Chapter II. The three tableaus of Figure 5.2.2 show the 


procedure of the equivalent algorithm. 


ale 











Figure 5.2.1 
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Figure 5.2.2 
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5.3 Revised Concept of a Primal AILP Algorithm 


Although the equivalent algorithm developed in the previous 
section is interesting from a theoretic standpoint, it is useless as 
a computational technique. It produces exactly the same results as 
the G-Y algorithm while requiring twice as much computation. Further- 
more, the equivalent algorithm would suffer from the same ineffi- 
ciencies, such as redundant cuts, that are inherent in the G-Y al- 
gorithm. The one redeeming factor in this equivalent algorithm is 
that it embodies the concepts essential to the handling of a non- 
integer tableau in a primal approach. Discussion of how these 
concepts may be employed to overcome some of the inefficiencies of 
the G-Y algorithm is reserved for the next section of this chapter. 
However, before pursuing this issue further, the present concept of 
a primal AILP algorithm must be revised. 

Recall from Chapter I that the primary need for a primal AILP 
algorithm was due to the generally unreliable character of the al- 
gorithms that were available. Computational experience with these 
algorithms showed that it was impossible to predict an upper limit 
for the number of iterations required for the solution of any parti- 
cular problem. This characteristic applies equally to the G-Y 
algorithm. However, the main virtue of this algorithm is that it 
produces intermediate feasible integer solutions during the computa- 
tion process. When computation costs prohibit a continuation to the 


optimum solution, the best of the intermediate solutions might be 
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useful. Indeed, a relatively ‘sooda' solution is better than no 
solution. Qn the other hand it is not strictly necessary for all 
intermediate solutions to be integer. The main concern 1s that 
feasible integer solutions should recur in the solution process and 
that the frequency of their recurrence can be predicted. It is in 
this context that the modified algorithm presented below will be 
referred to as a ‘primal’ method. All solutions are primal feasible, 


but not all of them are integer. 


5.4 Elements of a Modified Primal AILP Algorithm 


The AILP algorithm that is developed here is very similar to the 
equivalent algorithm developed in Section 5.2. Thus the rudimentary 
concepts of the different types of cuts will not be repeated. How- 
ever, several additional details must be supplied before the modifi- 
cations can be proposed. The first of these details cstablishes the 
capability of restoring a non-integer tableau to all-integer entries 
at any stage of the Simplex method: 

“If the initial tableau is all-integer, a Simplex solution 


procedure can be stopped at any stage of the computations and a 


series of Gomory cuts (5.1.2) can be employed to restore integrality. 


ayes 





Proof: 


Suppose that at some stage of the Simplex computations a matrix B 


has been inverted. The inverse of B can be expressed as: 


where: 


B = the adjoint of B (which must be integer) 


det B = determinant of B (product of pivot elements 1 through i-1) 


h 


Let p' = det B (in the it tableau) 


lo'| - absolute value of D' 
Every entry in tableau i can be expressed as: 


a; 5; = 4/0 H,D integer 


A cut (5.1.2) generated from tableau i would have a pivot element of 


the form: 
“a j = -(h/ |p" ) h,D integer 
0< h«< Ip) 


As a result of pivoting on -d. k in tableau i, tableau 7+ 17 would 


have 


Peery DU yep. Seen 
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Obviously since h is strictly less than p! by an integer amount, 
at most n & 0" | - | pivot steps could be required to begin with 
tableau i and produce an all-integer <ableau (because pitn = sag 

The above discussion seems to imply that one could start a Simplex 
solution of the LP problem and, by keeping track of D, decide to in- 
terdict with the cuts to produce an integer solution at an “opportune” 
moment. This is true, but unfortunately the integer solution could 
not be guaranteed to be feasible. Such an approach, though appeal- 
ing, 1S perhaps too ambitious. However, further restrictions can be 
applied to preserve feasibility. 

Suppose that the Simplex method is begun with an all-integer 
tableau. The Simplex procedure is allowed to continue until a non- 
integer solution results (integrality in the rest of the tableau is 
of no concern). Geometrically speaking this transition from integer 
vertex to non-integer vertex 1s pictured as a move along an edge of 
the convex polyhedron. Recall that movement along this edge has 
resulted from pivoting in a- certain tableau column k. Let ap be the 
integer column vector of values of the problem variables in tableau i 


and an’ the non-integer vector that results from pivoting in column k 


of tableau i. (ap 0 and an”) are the associated values of the objective 
function). 
"If a Simplex pivot in column k causes a transition from a pre- 


viously acquired integer solution in tableau i to a non-integer 


Solution in tableau i+], at most n = lo'* | - | pivots in column k 





on cuts (5.1.2) are required.to regain a feasible integer solution 


having the property das 2 a 


1 it 
9 


0° 
Proof: 
Consider the Simplex pivot from tableau i.to tableau i+] 
it] oi 1 
ag = 8g - > & 
where A is a positive scalar and a 1S a lexico-negative pivot 


column by virtue of the generalized Simplex pivot rules. In tableau 
i+] alt becomes lexico-positive. Successive pivots in column k 
on cuts (5.1.2) will change the.solution vector a, in the following 
way: 


i+j _ itj-1 _ i+j-1 
a ay ae oe (5.4.1) 


where a is a positive scalar and a, is always a-lexico-positive 


column vector. Thus (5.4.1) shows that the successive a, form a 
series of lexico-decreasing solution-vectors beginning with . 


i+m+] 


However, the values that ay may take are bounded from below (in 


<a) 


a lexico-graphic sense) by ay because a.,, + a.45 +..-+ Gi s 


i+] it2 
(the Gomory cuts are specifically guaranteed not to exclude the integer 


Jonnell i+m+] 


solution ap). This establishes that 4.0 = 49.9 and a = integer 


1 
m+1 <n and proves the proposition. 

Notice that the proof just given immediately allows one to avoid 
the redundant cuts encountered in the G-Y algorithm. When a cut in 


the G-Y algorithm would be generated from a constraint S. such that 


aya 





ay o/ yk = integer, the cut could be neglected and the constraint S, 
could be chosen as the pivot row. The next solution would be integer 
and the rest of the tableau non-integer. When later G-Y pivot steps 
produced a non-integer solution one could revert to cuts (5.1.2) to 
regain an integer solution (they can never be redundant). This is 
one concept of a modified algorithm. Such an algorithm would employ 
both of the types of cuts that have been discussed. In the general 
case the cuts (5.1.2) would no longer be equivalent to the cuts 
(5.2.2) as before. This is so because one would no longer generate 
them in the very specialized manner prescribed by the equivalent 
algorithm. However, there is another factor present that has led the 
uuthor t- conclude that the cuts (5.2.2) should be abandoned completely 
in favor of a generalization of the equivalent algorithm. This con- 


clusion 1s based on the fact that there is an alternative to the cut 


(5.1.2) that is more effective. 


5.5 Alternative Cuts in the Modified Algorithm 


Recall the geometric interpretation of the function of the cut 
(5.1.2). In any given tableau the current solution lies at the origin 
of the non-basic (parametric) coordinate system. The cut (5.1.2), 
when appended to that tableau, makes the origin of non-basic space 
infeasible (i.e., 1t cuts away the origin). The effectiveness of any 
cut depends largely on how much of the convex polyhedron it excludes. 


Therefore, one would like to generate a cut whose intersections with 


BAO. 





the various coordinate axes lie as far away from the origin as 
possible. One cut 1s uniformly ‘deeper’ than another cut when its 
intersections with the coordinate axes are uniformly farther from 
the origin. 

Hadley [1] has presented a cut, derived from (5.1.2), that is 
at least as deep and sometimes can be uniformly deeper than (5.1.2). 


He considered the cut: 


n 


ris adigt Do dy s)(-ty ae 
jl 


and showed that if there are any as having the property qi > qi 0? 
the term (ds -JAeds /(-ds 9) can be substituted for di in 
(5.5.1) to yield an alternative cut that is deeper in the t. direction 
for every column j in which the conditions are met. In the limited 
number of problems on which the modified algorithm has been tested, 
these deeper cuts have always accelerated the convergence process. 
However, these deeper cuts should be used with care. Modification of 
the cut (5.5.1) yields a constraint for which the slack r‘ is no 
longer restricted to integer values. To date the author has used a 
heuristic device to overcome this difficulty. Recall] that the term 

“ha can be expressed as ho/ p' , O< hp < lo"). Once the modi fica- 
tion has been applied to (5.5.1) the resulting cut has been multiplied 
by the scalar ({p"] = ho) to yield deeper cuts that have ‘worked’ on 


the problems tested. This is not to say that such a procedure wil! 
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restrict r' to integer values. Indeed, no theoretical basis has been 
established for this technique; it was developed strictly through 
trial and error. Actually the whole topic may be of little conse- 
quence because the inclusion of the modified cut simply converts an 
AILP problem into a MILP problem. Chapter VI presents an algorithm 
very similar to the modified algorithm that will handle the MILP prob- 
lem. The decision on how to handle these deeper cuts can only come 


with more experience than is presently available. 


5.6 Pivot Selection Rules for the Modified Algorithm 


There is one final point that should precede a formal statement 
of the modified algorithm. This point relates to a modification of 
the Simplex pivot rules that should be used to accompany this algorithm. 
Close scrutiny of the proofs in Section 5.4 will show that it is pos- 
Sible to make a complete cycle from an integer solution to a non-integer 
solution and back to the same integer solution. This process is called 
a minor cycle. A minor cycle of the modified algorithm is roughly 
equivalent to a zero cut in the G-Y algorithm. The possibility that 
this phenomenon can occur in the modified algorithm presents a problem 
for convergence that has not been overcome to date. However, computa- 
tional results on the Jimited number of problems that have been 
attempted have been very favorable in comparison to other AILP al- 


gorithms and no convergence difficulties have been experienced. 
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Surely one would want to avoid minor cycles whenever possible. 
It will be shown below that the impending occurrence of a major 
cycle can be recognized in the algorithm. This fact will motivate 
a change in the Simplex pivot selection rules so that major cycles 


are produced whenever possible. 


"If in tableau 71 a Simplex pivot is selected on row 7 in 


column k such that a; 0/4; KR 2 lo'l, a major cycle is eminent." 
> > 


Proof: 

Recall that every entry in tableau 1 can be expressed as H/ p! P 
_ 1 

H,D integer. Let a; 9/a; , =A, A> ‘Dale 


Thus: 


: i+ as 
In tableau i+m, ay We Ti = nase <n are lexico-positive. 


and successive cuts (5.1.2) yield a series of lexico-decreasing 
solution vectors ane by virtue of the expressions developed in 
Section 5.4. 


And: ° c ° 
: i¢m _ _i+] i+] 
y= 89 (O54 + O42 tee% Sam-1) 20 


However, notice that that is another feasible integer solution: 


a) = ay - 6[D'l a (5.6.1) 
6 = [a/|p'l] 


Hee 





that is lexico-larger tnan ay. From this we can deduce that 


(a. )<) -§|p"| 


+ : Tee ; 
1+] C542 = 


1+m 


Since the Gomory cuts are guaranteed not to exclude the integer solu- 
tion (5.6.1). Integrality must be restored inn < Ip" pivot steps 
and (5.6.1) must be the resulting integer solution. 

In order to give preference to Simplex pivot steps that are 
known to initiate major cycles one should select a lexico-negative 
column k for which ai 9/44 k > lo" if one is available. This is 
equivalent to the technique used in the G-Y algorithm to avoid zero 


cuts. 


5./ Computational Results 


Obviously the derivations of the previous sections show that 
there are several alternatives available in both cut generation and 
pivot selection. The following is a formal statement of the combina- 
tion of rules that have proved to be the most effective. (An all 
integer initial tableau is assumed. ): 

A. In tableau i select the most negative column k 

that has a most binding row 7 with mo > "|. a 0 
(if there are no negative columns, tableau i js 
optimal.) (If there is no as > 0 the problem is un- 


bounded. ) 





B, Perform a Simplex pivot step on the element dik? 
C. If tableau 1+] has a non-integer solution, select 
a particular row that displays a non-integer value 
and generate the cut (5.1.2). GO 1T00D. If tableau 
i+] has an integer solution, GO TO A. 

D. If the conditions described in Section 5.5 are met, 

derive a deeper cut. GO TOE. 

E. Pivot on the cut in column k. Erase the cut. G0 

OG. 

Table 5.6.1] 1S a summary of results of the modified algorithm 
when compared to other AILP algorithms. This table was produced in 
part from data reported by Gonzales [2]. The test problems are ones 
that are designed to be relatively hard [3]. 

An interesting fact that was noted in the test problems was that 
the cases where the modified algorithm shows a markedly better 9e;- 
formance over that of the G-Y algorithm was not limited to the prob- 
lems on which G-Y generates redundant cuts. This is primarily due to 
the fact that the deeper cuts discussed above sometimes prove to be 
much more efficient than the specialized cuts of the G-Y algorithm. 
However, it should be pointed out that the problems on which the 
algorithms were compared were designed (by Thompson [3]) to be 
difficult for the algorithms that were in existence at that time. 
Thus the tests cannot be considered to provide conclusive evidence 


in favor of any particular algorithm of the group. In view of the 


Peon 





TABLE 5.6.1 





| DIMENSIONS ITERATIONS REQUIRED | 
OF TABLEAU | ALGORITHM II | THOMPSON HOPPER 


6 
25 
3 


* Indicates that the program was stopped after this number of iterations 


without having found the answer. 


All of the algorithms require roughly the same amount of time to perform 


an iteration. They all consist of Simplex pivot steps. 
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fact that certain of the available algorithms give their best 
performances on specific classes of problems [4], many more 

tests on practical problems would have to be performed before 
judgements can be made on the suitability of the modified algorithm. 


The author's opinions on this topic are reserved for Chapter VII. 





— a 4 7 
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CHAPTER VI 


ESSENTIALS OF A PRIMAL MILP ALGORITHM 


This chapter takes the concepts behind the modified AILP as a 
base and generalizes them only slightly in order to outline the 
essentials of a primal MILP algorithm. The cuts that are generated 
in this algorithm were developed first by Gomory in his Algorithm III 
for the MILP problem (a dual feasible method). The approach here is 
exactly the same as was taken in developing the modified AILP algorithm. 
Therefore, the geometric interpretations and background materia! are 
dispensed with here in favor of a simple statement of the algorithm. 
An example problem is presented that takes the MILP through several 
Steps of its procedure. Little computational experience with this 


algorithm is available and none is reported. 


6.1 Relationships Among LP, AILP and MILP Problems. 
Recall from Chapter I that both types of the ILP problem are forn- 


ulated in the same way as an LP problem; the only difference being the 
case of ILP problems. Further, it should be obvious that al! cutting 
plane algorithms approach ILP problems in the same manner. They all 
attempt to add new constraints to reduce the solution space. There 
are different types of Simplex algorithms (e.g., primal Simplex, dual 


Simplex), but every cutting plane method is, in some way, a modification 
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to one of these Simplex algorithms. They all solve the ILP problem 
by modifying the associated LP problem through the use of cuts. 

In Chapter V it was seen that cuts from two seemingly very 
different algorithms (G-Y and Algorithm I) were actually rather 
strongly related. In the final statement of the modified AILP algo- 
rithm, cuts from Algorithm I were employed in a primal context and 
the specialized cuts of the G-Y algorithm were ignored completely. 
When one realizes that Algorithm I for the AILP problem can be con- 
sidered to be a special case of Algorithm III for the MILP problem 
(the subset of integer-restricted variables contains the entire set 
of variables), the possibility of converting Algorithm iII to a 
primal approach becomes evident. In fact the approach remains the 
same, the only difference lies in the way the cuts are generated, and 


not the way they are used. 


6.2 Cuts for the Primal AILP Algorithm 


The cuts that are generated in Algorithm III represent a generali- 
zation of the cuts generated in Algorithm I. The two algorithms are 
completely parallel except for the cut generation. Therefore, Algorithm III 


will not be reviewed; its cut generation procedure wil! be presented. 


Consider a row of the tableau corresponding to some integer-constrained 


problem variable 


n 
BaP Ga + » a, y(-t,) (oecel) 





Let: 


R be a set composed of all non-basic (parametric) variables ue j=li sean 
R,ER contains ty for which a; , > 0 and t, not integer-cons trained 


R €R contains ty for which a.» <0 and t, not integer-constrained 


k 
I €& contains ty for which ty is integer-constrained 


Gomory considered a tableau row (6.2.1) for which Xi. is constrained 
to integer values and showed that feasible solutions that satisfy (6.2.1) 


also satisfy: 


n 


r' = “9. ,0 if 3 (-9, sje? (6.2.2) 


° 


j=l 


where the 9. ; are selected to produce the deepest cut by the following 


rules: 
“53 049 + 
d. ,0 qa. aj /(1-d, 0) > jER_ 
“a d j€l and d. .<d 
ee od tw 


ds gll-d; j)/(-d; 9) » jl andd; , > d; 4 


The cut (6.2.2), though designed for a dual feasible method, can be used 
at any stage of a Simplex method in the same way as cuts from Algorithm I 


were employed in Chapter V. However, it should be noted that r’ in (6.2.2) 
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iS not an integer-constrained variable. If the parametric variable ty 
is a non-basic proxy for an integer-constrained problem variable; t 
is also integer-constrained. If a Simplex pivot step is performed in 
tableau i in column k on a constraint S. and the slack S. 1s not 
integer-constrained, then ty in tableau i+] will become a non-basic 
proxy for Ss; and as such wil! no longer be integer-constrained. Thus 
it will be necessary to keep track of the Simplex pivot steps in order 
to decide to which subset of R (R,5 R_, I) a particular parametric 


variable belongs. 


6.3 Statement of the Primal MILP Algorithm 


Attention is recalled to the geometric interpretation of the Sim- 
plex pivot steps. This interpretation will not be repeated here. 
However, notice that the concept of stopping a Simplex method at some 
stage and cutting back to a (ILP) feasible solution is just as valid 
in the MILP case as it was seen to be in the AILP case. Therefore, 
the statement of the primal MILP algorithm exactly parallels the state- 
ment of the AILP algorithm. (Assume a starting basic feasible solution 


to the MILP problem. ) 


A, Select a column k for pivoting by the normal Simplex 
pivot selection criteria. Pivot in column k. GO JO B. 

B. If the resulting tableau produces values of the 
problem variables that satisfy the integrality con- 


straints, GO T0 A. Otherwise, GO T0 C. 





C. Select some row of the tableau that represents 
an integer-constrained variable X where Ks does 
not have an integer value in this tableau. From 
this row generate the cut (6.2.2). Pivot on this 


cut in column k. Erase the cut. GO TO B. 


This algorithm obviously is intended to terminate in Step A when the 
normal Simplex pivot selection criteria indicate either an optimal 
or an unbounded solution. 

Though the procedure described above parallels the procedure of 
the modified AILP algorithm presented in Chapter V, it has one some- 
what weaker point. When a Simplex pivot step transitions from a 
previously attained (MILP) feasible solution to one that is not (MILP) 
Fee ibie, no precise statement can be made about the maximum number 
of cuts that will be needed to regain a (MILP) feasible solution. 
However, it is known that this will be accomplished in a finite number 
of iterations because Gomory proved a related supposition in his proof 
of convergence for Algorithm III]. Without resorting to his detailed 
proof, the proposition may be supported by noting that column k (the 
pivot column) becomes lexico-positive as a result of the Simplex pivot 
Step that causes the transition. Every successive pivot step that is 
performed in column k on a cut (6.2.2) produces a Jexico-smaller solu- 
tion and leaves onan k lexico-positive. Thus the series of cuts 
produces a series of lexico-decreasing solutions that is bounded from 


below by the previous (MILP) feasible solution (since the cuts 





specifically cannot exclude this solution). Gomory showed that a 
Sequence of lexico-decreasing solutions similar to the one described 


here will ultimately produce the desired (MILP) feasible solution. 


6.4 An Example Problem 


Since this chapter has provided only a sketchy outline of the 
essentials of the primal MILP algorithm, and since the author's 
experience with this technique is very limited, the example problem 
presented below will illustrate the details of this algorithm. The 
problem that is under consideration is the same as the one used 
throughout this paper to demonstrate the AILP algorithms (though the 
objective function has been changed). However, the variable Xo is 


now allowed to take on any non-negative value: 


max OX + 2X5 


A 
on 


~Xy + oX5 < 
xX, + XS 9 
Xi» Xp > Oars xy integer 
Figure 6.4.1 shows the initial tableau and the tableau resulting from 
a normal Simplex pivot in column 2. (The most negative column was not 


chosen in this caSe because this allows the MILP algorithm to be demon- 


strated.) Notice that Xx, has an integer value (0) so that this solution 


acy 





is (MILP) feasible. The first tableau of Figure 6.4 2 is the result 
of another normal Simplex pivot step. However, notice that the in- 
teger restrictions on X, are violated. Therefore, a cut (6.2.2) is 


generated from the constraint row corresponding to x Pivoting on 


1° 
this cut restores (MILP) feasibility and in this case with a better 
value of the objective function than was previously available. The 
fact that Xo also assumes an integer value has no bearing whatsoever 


on the algorithm, since it is a=lowed to take on any non-negative 


value. 












nee 


Figure 6.4.2 








CHAPTER VII 


SUMMARY, CONCLUSIONS, AND SUGGESTED EXTENSIONS 


7.. 1 Summary 


The work reported here was composed of two distinct, though 
related, parts. The first of these was confined to a presentation 
of the G-Y algorithm with an accompanying discussion of the manner 
in which this algorithm has been coded into a problem-oriented lan- 
guage. Some details of the data structure and data management poli- 
cies employed in this POL were discussed. Complete program docu- 
mentation was included as an appendix to this report. Finally, 
detailed instructions on the use of this POL were given through com- 
plete documentation of the commands and their uses, an illustrative 
example problem encompassing all of the OPAILP’s capabilities, and 
several random tips on some practices that should be avoided while 
using the system. 

The second part of this paper became involved in the topic of 
primal ILP algorithms in general. Some geometric interpretations 
accompanied an explanation of a specific problem that leads the G-Y 
algorithm to experience unnecessarily slow convergence rates, namely 
the generation of redundant cuts. The first step in an attempt to 
show that these problems could be overcome was an illustration of 


how cuts taken from Gomory’s Algorithm I could be employed in a primal 





feasible context to produce an AILP algorithm that was equivalent to 
the G-Y algorithm. This equivalent algorithm demonstrated a capability 
of working with a completely non-integer tableau while producing 

primal feasible integer solutions to the AILP problem. The concept 

of a primal feasible AILP algorithm was altered slightly to include 

the case where integer solutions could be guaranteed to recur at vari- 
ous intervals in a modified Simplex algorithm. The cuts of Algorithm I 
were employed in a special! way to produce a modified primal AILP algo- 
rithm. It was established that this algorithm could be guaranteed to 
produce primal feasible integer solutions at well defined intervals in 
the solution process. The performance of this algorithm was compared 
to various other AILP algorithms on a limited number of test problems. 
The results of these tests were deemed to be inconclusive. 

Finally, the concepts of the modified AILP algorithm were ex- 
tended to produce an outline of the essentials for the first primal 
feasible algorithm for the MILP case. The procedures of this algorithm 
were illustrated through an example problem, though no other computa- 


tional results were offered. 


7.2 Conclusions 

No attempt has been made to compare the computational efficiency 
of the computer programs developed here to other existing packaged 
programs for the solution of AILP problems. The reason for this is 


that no other programs based on the G-Y algorithm are publically 
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available. As a result, comparisons among the various programs would 
be highly algorithm-dependent and could not be construed to measure 
the relative efficiency of the coding involved. 

Along these same !/ines the author has specifically avoided 
giving estimates of the number of G-Y iterations that can be expected 
for problems of various sizes. While such estimates may be properly 
offered as rough guidelines in terms of LP problems, no such estimates 
are valid in terms of AILP problems. Factors other than problem size 
are far more relevant to the prediction of the number of iterations to 
be expected, namely the relative primeness of the coefficients in the 
various constraint rows. 

As a result of the above-mentioned considerations, there remains 
no basis for judgements about the relative efficiency of the coding 
employed to develop the POL reported in the first part of this paper. 
The timing tests and solution attempts on a wide variety of problems 
that would provide the information needed to make these judgements 
have been subjugated to the author's more active interest in the 
work reported in the second part of this paper. Time constraints have 
not permitted the two topics to be pursued to the same extent and the 
former has suffered as a result. Therefore, OPAILP should be con- 
Sidered an experimental component of OPTECH until further computationa! 
experience can be gained. 

It was reported tn Chapter V that performance tests with the modi- 


fied algorithm were too limited to be conclusive. On some of these 





problems the modified algorithm was impressively more efficient than 
all other algorithms tested. However, it should be pointed out that 
the performance of each of the existing algorithms has always proved 
to be strongly influenced by the type of AILP problem that was being 
tested (i.e., some algorithms are very efficient on cutting stock 
problems and inefficient on covering problems, etc.). The author 
believes that this characteristic will be common to the modified 
algorithm. There are undoubtedly many problems that will yield more 
readily to the G-Y approach than to the modified algorithm. This 1s 
so because the modified algorithm was designed to cope with exactly 
those problems that are troublesome for the G-Y algorithm. For the 
problems on which G-Y is efficient, the modifted algorithm wil] 
probably be relatively inefficient because cuts cannot be generated 
in the modified algorithm more often than in every other tableau. 

If the G-Y cuts on a particular problem happen to be very efficient 
the modified algorithm will suffer by comparison. In any event, many 
more tests are required before the modified algorithm can be placed 
in its proper prospective with regard to the other AILP algorithms. 
In the absence of supporting computational evidence, the author 
believes that the modified algorithm derives some importance from its 
instructional value. Its development bridges the gap between the G-Y 
primal algorithm and Gomory’s dual feasible Algorithm I and shows that 
the seemingly different cuts employed in these two algorithms are, in 


fact, very strongly related and in some case may be equivalent. 
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The outline of the essentials of a primal MILP algorithm that 
was presented here may seem to be only an obvious extension of the 
modified algorithm. In fact this 1s true, but the author feels that 
this may prove to be by far the most important topic discussed here. 
MILP formulations abound in all areas of managenient and engineering, 
but the complexity of these problems is often directly proporticna|! 
to their economic and social importance. Too often the proper forn- 
ulation of these problems has been abandoned in favor of an approach 
that is computationally more feasible, though less precise. There- 
fore, the importance of developing a computational feasible MILP 
algorithm cannot be overemphasized. If that algorithm could be 
primal feasible as well, it would be all the more valuable. The 
author does not claim that the MILP algorithm outlined here represents 
the answer to this problem. However, the need for an efficient MILP 


algorithm is too pressing to allow any possibility to go untested. 


7.3 Suggested Extensions 


Work on every area encompassed by the current research has left 
several questions unanswered. These questions point to the obvious 
need for further work. However, the author feels that continuation of 
this work can be organized in a particular way to provide maximum re- 
turn. 

First the OPAILP component of OPTECH should be thoroughly tested 


in order to insure that the coding is precisely the best for the 
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Situation. This testing should be performed on a variety of types 
of AILP probelms; specifically problems that have already been used 
to test other AILP algorithms. This will allow for a determination 
of the types of problems on which G-Y is to be preferred over some 
other algorithm. 

Secondly an attempt should be made to interface the LP routines 
and the AILP routines of OPTECH to allow for the interaction necessary 
to solve a MILP problem via Benders’ algorithm. OPAILP has been 
incorporated into OPTECH in such a way as to allow the requirements 
for data transfer through the subsystem's COMMON area. 

The need for computational experience with both the modified AILP 
algorithm and the MILP algorithm is obvious. Investigations into the 
proper pivot selection and cut generation strategies for these algo- 
rithms should be more extensive. Testing of the MILP algorithm on 
(0,1) decision variable problems is specifically recommended. Finally, 
both of these algorithms lack a proof of finite convergence. Serious 
attempts to provide these proofs would more than likely provide the 


proper pivot selection and cut generation strategies as a by-product. 
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APPENDIX A 
PROGRAM DOCUMENTATION 


This Appendix contains ICETRAN listings of all of the programs 
that make up the CDB's for OPAILP. Each listing accompanied by an 
explanation of what the program does and how it fits into the se- 
quencing of the algorithm. These explanations are contained on the 
side of the paper that faces the listing itself to facilitate an 


easy view of the program that is being discussed. 
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Documentation of Subroutine OPMANX 
OPMANX does the following things: 


1. Sets up alphanumeric data to guide the input of the 
tableau,L,G,S000, » OBJ, etc. All stored in LES 


2. Computes the storage space required for the problem 
and stores this in LES(6) 


3. Defines the approiate dynamic arrays to the sizes 
called for by data from the AILP command. NCL=no. of 
cols NRO=no. of rows of the problem. 


4, Finally sets a flag that will let the RIGHT SIDE 
command be processed. 
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ICETRAN LISTING OF SUBROUTINE OPMANX 


20 


SUBROUTINE OPMANX 

COMMON IAM(P),B(P),NAMR(P),A(P),X(P),NAMC(P), IROF(P) 
COMMON Y(P),E(P),P(P),ICLF(P),KM(50) ,KBCD(9) ,KP(10) 
COMMON Z(15),1(6),KIT, IRW,LR,LK,LC,KPIV,KZZ,1IS,IPZ,KZ 
COMMON NCAS(10) ,KL,IE,NCL,NRO,LES(10), IPIV,LEC 
DYNAMIC ARRAY IAM,NAMR,NAMC, IROF,ICLF 
LES(1)=1087816401 

LES(2)=1077952723 

LES(3)=1077952711 

LES(4)=1077952576 

LES(5)=-487526160 

LES(6)=4*( (NCL*(NCL+NRO+6) )+(3*NRO)+7) 

DEFINE IROF,NCL+NRO+1 ,STEP=5 

DEFINE NAMR,NCL+NRO+1 ,STEP=5 

DEFINE NAMC,NCL+1 ,STEP=1 

DEFINE ICLF,NCL+1 ,STEP=1 

DEFINE IAM,NCL+1,POINTER,STEP=1 

IPT=NCL+1 

DO 20 JUX=1,IPT 

ICLF(JJX)=0 

DEFINE IAM(JJX),NCL+NRO+1 ,STEP=5 

IPZ=-97 

RETURN 

END 
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Documentation of subroutine OPRHSI 


OPRHSI is entered when the RIGHT SIDE command is encountered. 
It does the following things: 


1.Checks for proper command sequencing. 

2.Reaas a sequence of cards noting the type of constraint 
(in IROF),the name of the constraint (in NAMR), and stores 
the coefficient value in the first column of the tableau(IAM). 


3.The row names are used later in OPTAB to guide in placing 
the matrix coefficients in the proper place. 
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ICETRAN LISTING OF SUBROUTINE OPRHSI 


2 | 


ND 


24 
23 


oy on 


o™! 


100 
10] 
950 
102 


SUBROUTINE OPRHSI 

COMMON IAM(P),B(P),NAMR(P),A(P),X(P),NAMC(P), IROF(P) 
COMMON Y(P),E(P),P(P),ICLF(P) ,KM(50) ,KBCD(9),KP(10) 
COMMON Z(15),T(6) KIT, IRW,LR,LK,LC,KPIV,KZZ,IS,IPZ,KZ 
COMMON NCAS(10),KL,IE,NCL,NRO,LES(10),IPIV,LEC 
DYNAMIC ARRAY IAM,NAMR,NAMC, IROF,ICLF 
IF(PZ+97)21,22,21 

WRITE(6,950) 

WRITE(6,102) 

CALI EXL | 

IRW=2 

READ(5,100) ITRM,ISLK,NRT,IAT 
IF(ITRM-LES(4))2,3,2 

IAM(1,1)=0 

IROF(1)=0 

IRW=IRW- 1 

NAMR(1)=LES(1) 

IPZ=-99 

RETURN 

IR 72+3)23,24.23 

WRITE(6,100) ITRM,ISLK,NRT, IAT 

NAMR(IRW)=NRT 

TAM(1,IRW)=IAT 

IFCISLK-!FS(2))5,4,5 

IROF(IRW)=1 

GO TO 1U 

IF(ISLK-LES(3))7,6,7 

IROF(IRW)=0 

TAM(1, IRW)=IAM(1,IRW)*(-1) 

GO TO 10 

IF(ISLK-LES(4))9,8,9 

IROF(IRW)=-1 

IRW=IRW+1 

GO TO 1 

WRITE(6, 100) ITRM,ISLK,NRT,IAT 

WRITE‘6,101) ISLK 

WRITE(6,102) 

CALL EXIT 

FORMAT (3A4,2X ,16) 

FORMAT(//,6X,'IN ABOVE CARD-' ,A4,'IS ILLEGAL. ‘) 
FORMAT(//,6X,°COMMANDS OUT OF SEQUENCE. ') 
FORMAT(//,6X,'JOB FLUSHED’ ) 

END 
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Documentation of subroutine OPTAB 


Control! enters OPTAB when the TABLEAU command is encountered. 
It does the following things: 


1. It reads a sequence of data cards terminated by a non-blank 
field 1-4. 


2. It searches the list of row names that was input in the 
OPRHSI subroutine. It assumes that the coefficients of any 
column are input together. 

3. It names the columns according to the user supplied names. 


4, It changes greater than or equal to inequalities to less 
than or equa! to. 


5. Once the sequence has ended, OPTAB fills out the rest of 
the tableau with the -I1 matrix representing the problem variables. 


6. Various constants are set, and contro! is returned to the 
command processor. 
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ICETRAN LISTING OF SUBROUTINE OPTAB 


40 


4| 


nS —_ 


58 
59 


2] 


12 
3 


SUBROUTINE OPTAB 

COMMON IAM(P),B(P),NAMR(P),A(P),X(P) ,NAMC(P) , IROF(P) 
COMMON Y(P),E(P),P(P),ICLF(P),KM(50) ,KBCD(9) ,KP(10) 
COMMON Z(15),1T(6), KIT, IRW,LR,LK,LC,KPIV,KZZ,1S,IPZ,KZ 
COMMON NCAS(10),KL,TE,NCL,NRO,LES(10) ,IPIV,LEC 
DYNAMIC ARRAY IAM,NAMR,NAMC,TROF,ICLF 
IF(IPZ+99) 40,41 ,40 

WRITE (6,950) 

WRITE(6, 102) 

CALL EXIT 

KT=0 

WRITE (6,105) 

KIT=2 

READ(5,100) ITRM,NCCK,NRCK, ITEL 
IF(ITRM-LES(4))2,3,2 

JAX=IRW+KIT=1 

TS=IRW+] 

DO 21 MO=1,KIT 

DO 58 KOP=IS,JAX 

TAM(MO,KOP )=0 

IF(MO-1)59,21,59 

MI=MO+IRW- 1 

NAMR (MI) =NAMC (MO) 

IROF (MI )=0 

TAM(MO, MI )=-1 

RELEASE IAM(MO),LOW 

TE=MI 

IPZ=0 

KPT V=0 

NAMC (1)=LES(1) 

TRW=MI 

RELEASE NAMR ,LOW 

RELEASE NAMC,LOW 

LEC=49 

RETURN 

IF(LR#2)42,43,42 

WRITE(6,100) ITRM,NCCK,NRCK,ITEL 
IF(NCCK)5,4,5 

WRITE(6,100) ITRM,NCCK,NRCK,ITEL 
WRITE(6,101) 

WRITE(6,102) 

CALL EXIT 

IF(KT)6,10,6 

IF(NAMC(KIT)-NCCK)7,12,7 
IF(NRCK)14,13,14 
WRITE(6, 100) ITRM,NCCK,NRCK, ITEL 
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ICETRAN LISTING CONT INUED-OPTAB 


14 


V7 
19 
18 


wd 


100 
101 
102 
103 
104 
105 
950 


WRITE (6,103) 

WRITE (6, 102) 

CALL EXIT 

DO 15 JXX=1,IRW 

IF (NAMR(JXX)-NRCK)15,17,15 
IF(1ROF(JXX))18,19,18 

ITEL=ITEL*(-1) 

TAM(KIT,JXX)=ITEL 

GO TO 1 

CONTINUE 

WRITE(6,100) ITRM,NCCK,NRCK,ITEL 

WRITE(6,104) NRCK 

WRITE(6,102) 

CAEL | 

RELEASE IAM(KIT),LOW 

KIT=KIT#1 

NAMC (KIT )=NCCK 

KT] 

GO TO 14 

FORMAT (3A4,2X,16 

FORMAT(//,6X,'IN ABOVE CARD, COL. NAME IS REQ.') 
FORMAT(//,6X,'JOB FLUSHED' ) 

FORMAT(//,6X,'IN ABOVE CARD, ROW NAME IS REQ.') 
FORMAT(//,6X,'NO MATCH FOR-' ,A4,‘=-IN ABOVE CARD.') 
FORMAT (//) 

FORMAT(//,6X, COMMANDS OUT OF SEQUENCE‘) 

END 
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Documentation of subroutine OPCHOZ 


Control enters OPCHOZ when the ITERATE command is encountered. 
It does the following things: 


1, It checks to see if there are any equations in the formulation. 
2. If there are equations, OPCHOZ locates their row subscripts 
one at the time(LR). Control is passed to OPEQEL to have the 
equation taken care of. 

3. When there ore no more equations in the tableau, control is 


passed to OPPSLI where the solution process is continued until 
termination. 


-Ill- 





ICETRAN LISTING OF SUBROUTINE OPCHOZ 


40 


4] 


DN _ OO > Wm 


102 
103 
950 


SUBROUTINE OPCHOZ 

COMMON IAM(P),B(P),NAMR(P),A(P),X(P),NAMC(P),IROF(P) 
COMMON Y(P),E(P),P(P),ICLF(P),KM(50),KBCD(9) ,KP(10) 
COMMON Z(15),T(6),KIT,IRW,LR,LK,LC,KPIV,KZZ,IS,1PZ,KZ 
COMMON NCAS(10) ,KL,IE,NCL,NRO,LES(10),IPIV,LEC 
DYNAMIC ARRAY ITAM,NAMR,NAMC, IROF,ICLF 

IF(IPZ)40, 41,40 

WRITE(6,950) 

WRITE(6, 102) 

CALL EXIT 

KL=2 

KZ=KIT 

GO TO 4 

CALL OPEQEL 

IF(KL-IRW)9,9,10 

DO 2 NZ=KL,IRW 

IF(IROF(NZ))1,2,2 

LR=NZ 

GOO. 5 

CONTINUE 

RELEASE IROF,LOW 

WRITE(6,103) KPIV 

CALL OPPSL1 

RETURN 

FORMAT(//.6X,'JOB FLUSHED’ ) 

FORMAT(//,6X,'FEASIBLE ITERATION’ ,2X,14) 
FORMAT(//,6,'COMMANDS OUT OF SEQUENCE. ') 

END 
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Documentation of subroutine OPEQEL 


Control enters OPEQEL from OPCHOZ with a row subscript of 
an equation in LR. OPEQEL does the following: 


1. Finds the smallest coefficient in the row LR(absolute value). 
2, It puts the column subscript of the coefficient in LC. 


3. Calls the computation routine OPCMP1] to have a pivot performed 
on IAM(LC,LR). 


4, When the equation has been satisfied and there is but one 


coefficient in the row, this column is shifted to the last column 
and is no longer a part of the problem. 
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ICETRAN LISTING OF SUBROUTINE OPEQEL 


99 


— Ww OF MP 


So 


V7 


77 
16 
ie 


101 
102 
200 


COMMON 1AN(P) .B(P 

COMMON Y(P),E(P) 

COMMON 2(15),7(6 
) 
I 


) ,NAMR(P) ,A(P),X(P) ,NAMC(P) , IROF(P) 
,P(P), ICLF(P),KM(50) ,KBCD(9) ,KP(10) 

) KIT, IRW,LR,LK,LC,KPIV,KZZ,1S,IPZ,KZ 
KL,IE,NCL, NRO ,LES(10), IPIV, LEC 


E( 
COMMON NCAS(10), 
AM, NAMR,NAMC, IROF, ICLF 


DYNAMIC ARRAY 
KPOZ==1 
ILO=32000 

DO 1 JX=2,KIT 
IF(IAM(JX,LR))5, 
ITST=IAM(JX,LR)* 
GO T0 6 
ITST=IAM(JX,LR) 
TF(ITST- IL0)3 34 
ILO=ITST 

LC=JXx 
KPOZ=KPOZ+1 
RELEASE IAM(JX),LOW 

IF(KP0OZ)15,8,11 

IF(ITAM(1.LR))13,10,13 

IROF(LR)=1 

KL=LR+1 

SWITCH(IAM(LC), TAM(KIT) ) 

RELEASE IAM(LC),LOW 

RELEASE IAM(KIT),LOW 

KIT=KIT=1 

RETURN 

IPIT=IAM(1,LR)/IAM(LC,LR) 
IF(IAM(1,LR)-(IPIT*(IAM(LC,LR))))17,11,17 

NOR=NAMR(LR) 

WRITE(6,101) NOR 

WRITE(6, 102) 

CALL EXIT 

CALL OPCMP 1 

GO TO 99 

NOR=NAMR (LR) 

WRITE,6,200) NOR 

GO TO 77 

FORMAT(//,6X,'NO INTEGER SOLUTION EXISTS TO EQUATION’ ,A4,'.') 
FORMAT(//,6X,J0B FLUSHED’ ) 

FORMAT(//,6X, INVALID EQUATION, ALL COEFFICIENTS ARE ZERO.',/) 
END 


leg 
(-1) 


4 
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Documentation of subroutine OPPSL1 


Contro! enters OPPSL1 when there are no longer any equations 
in the tableau. There may or may not be another infeasibility. 
OPPSLI does the following things: 


1. It treats Phase I and Phase II exactly alike. It will find 

an infeasible row and make it a temporary objective function. 

It also considers for pivoting only those rows that are currently 
feasible. 


2, It will reject a zero cut as long as there are other negative 
columns to be examined. If it must take a zero cut, it will take 
a row for which aso # 0 where possible. 


3. It also keeps a subproblem of columns that have previously 
generated zero cuts. (see Gonzales convergence heuristics) 


4, Control continues to pass between OPPSL1 and OPCMP1] until 
the problem is finished through optimality or through an 
overrun of the iteration limit. 


5. The current objective function is always stored in ICLF 
in order to avoid rowwise referencing through the tableau. 
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ICETRAN LISTING OF SUBROUTINE OPPSL] 


102 


103 
101 
104 


200 


202 
201 


wd 


96 
og 


98 


89 


SUBROUTINE OPPSL1I 

COMMON IAM(P),B(P),NAMR(P),A(P),X(P),NAMC(P) ,IROF(P) 
COMMON Y(P),E(P),P(P),ICLF(P),KM(50,KBCD(9),KP( 10) 
COMMON Z(15),7(6) KIT, IRW,LR,LK,LC,KPIV,KZZ,1S,1IPZ,KZ 
COMMON NCAS(10),KL,IE,NCL,NRO,LES(10) ,IPIV,LEC 
DYNAMIC ARRAY IAM,NAMR,NAMC ,TROF,ICLF 
DIMENSION IOPT(1) 

EQUIVALENCE (QQDUB(1),IOPT(1)) 

LK=0 

LOOK=-] 

KZZ=] 

12=(1FIND( IAM, 1)+QQICOM)/4 

KZ=0 

KL=2 

LS=0 

LF=0 

LR=0 

IF(LK-1)3,4,4 

DO 6 JX=2,IRW 

T22=] 2+ JX 

IF(IOPT(122) )5,6,6 

LK=JX 

GO TO 4 

CONTINUE 

LK=1 

GO TO 200 

IF (LOOK-LK)200,201, 200 

LOOK=LK 

DO 202 JO=1,KIT 

ICLF(JO)=IAM(JO,LK) 

RELEASE IAM(JO),LOW 

DO 7 JX=KL,KIT 

KI=JX 

BECIGPrC IX) )8,7,7 

CONTINUE 

IF (KL=2) 46,96, 46 

IF(LK-1)99,98,99 

WRITE(6,111) 

WRITE (6,105) 

CALL EXIT 

WRITE(6,112) 

CALL OPAUT} 

RETURN 

11=(IFIND(1AM,KI)+#QQICOM)/4 
IF(LS)88,89,88 

LC=KI 
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ICETRAN LISTING CONTINUED-OPPSLI 


LS=] 

88 DO 12 JX=2,IRW 
IF(JX-LK)13,12,13 

13 I11=114+JX 
122=12+JX 
IF(IOPT(122))12,15,15 

15 IF(IOPT(111))12,12,16 

16 LF=LF+1 
1F(LF-1)9,10,9 

10 LO=IOPT(122)/IOPT(111) 
LR=JX 
GO TO 12 

9 ITST=10PT(122) /IOPT(111) 
IF(LO-ITST)12,17,18 

18 LO=ITST 
LR=JX 
GO TO 12 

17 ITT=I0PT(122)-(IOPT(111)*ITST) 
ITO=IAM(1,LR)-(IAM(KI,LR)*LO) 
IF(ITO=1TT)22,19,12 

IF(ITO+#ITT)20,12,20 

22 LO=ITST 
LR=JX 

20 IF(KZ)12,21,12 

21 LC=KI. 
KZ=] 

12 CONTINUE 
IF(LF)30,31, 30 

31 IF(LK=1)32, 35,32 

35 WRITE(6,110) 
WRITE(6,105 ) 
CALL EXIT 

32 LC=KI 
LR=LK 
CALL OPCMP1 
GO TO 102 

30 IF(LS=1)42,49,42 

42 IF(LC-KI)41,49,41 

49 LS=LR 

41 IF(LO)43,44,43 

44 IF(KI-KIT)45,46, 46 

46 LR=LS 
IF(LC-KZZ)47, 47,48 

48 KZZ=KZZ+1] 
SWITCH(IAM(LC) , IAM(KZZ) ) 
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ICETRAN LISTING CONTINUED-OPPSL1 


ITE=JCLF(LC) 
ICLF(LC )=ICLF(KZZ) 
ICLF(KZZ)=ITE 
LC=KZZ 
47 CALL OPCMP] 
GO 70 103 
45 KL=KI+] 
RELEASE IAM(KI),LOW 
GO TO 104 
43 LC=KI 
CALL OPCMP } 
IF(LK-1)50,1,50 
50 IF(IAM(1,LK))1,102,102 
105 FORMAT(//,°JOB FLUSHED' ) 
110 FORMAT(//,6X,'THIS PROBLEM IS UNBOUNDED. ' ) 
111 FORMAT(//,6X,'THIS PROBLEM HAS NO FEASIBLE SOLUTION') 
112 FORMAT(//,6X,’OPTIMALITY HAS BEEN PROVED. ') 
END 
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ICETRAN LISTING OF SUBROUTINE OPCMP1 


On 


CoO NSN WPS 


96 


15 
16 


— 


12 


— 
al 


99 


SUBROUTINE OPCMP1 

COMMON IAM(P),B(P) ,NAMR(P),A(P),X(P),NAMC(P), IROF(P) 
COMMON Y(P),#(P),P(P),ICLF(P),KM(50,KBCD(9) ,KP(10) 
COMMON Z(15),T(6),KIT,IRW,LR,LK,LC,KPIV,KZZ,1S ,IPZ,KZ 
COMMON NCAS(10),KL,IE,NCL,NRO,LES(10) ,IPIV,LEC 
DYNAMIC ARRAY IAM,NAMR,NAMC,IROF,ICLF 
DIMENSION IOPT(1) 

EQUIVALENCE (QQDUB(1),I0PT(1)) 
IF(IAM(LC,LR((5,4,4 

IDIV=-IAM(LC,LR) 

GO TO 98 

IDIV=IAM(LC,LR) 

DO 1 KX=1,KIT 

IF(KX-LC)2,1,2 

IF(IAM(KX,LR))3,15,7 
KPO=(IAM(KX,LR)-IDIV+1)/IDIV 

GO TO 8 

KPO=IAM(KX,LR)/IDIV 

IF(KPO)8,15,8 

CALL ALOCAT (IAM,KX) 

CALL ALOCAT (IAM,LC) 
12=(IFIND( IAM, KX )+QQICOM)/4 
I1=(IFIND(IAM,LC)+QQICOM)/4 

DO 9 JO=1,IRW 

122=12+J0 

111=11+Jd0 

1OPT(122)=IOPT(122)-( 10PT(111)*KPO) 
ICLF (KX )=ICLF(KX)-FCLF(LC)*KPO) 
IF(LES(6)-16000,3,1,96 

RELEASE IAM(KX),LOW 

GO TO 1 

IF(KX-1)6,16,6 

IPZ=IPZ+] 

GO TO 6 

CONTINUE 

IF(IAM(LC,LR))11,11,10 

DO 12 KX=1,IRW 

111=114+KX 

LOPT(I11)=IOPT(111)*(-1) 
ICLF(LC)=ICLF(LC)*(-1) 

KPIV=KPIV+1 

RELEASE IAM(LC),LOW 
IF(KPIV-IPIV)97,99,97 

WRITE(6,110) 

WRITE(6,1171) 

WRITE(6, 112) 
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ICETRAN LISTING CONTIUED-OPCMP 1] 


CALL OPAUT | 
CALL EXIT 
97 RETURN 
110 FORMAT(//,6X,'USER ITERATION LIMIT EXCEEDED. ') 
111 FORMAT(//,6X, CURRENT SOLUTION FOLLOWS. ') 
112 FORMAT(//,6X,'THEN JOB TERMINATED. ') 
END 
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Documentation of subroutine OF CMP] 


Control enters OPCMP1] when an updating operation is to be per- 
formed. The subscripts of the pivot element are LC,LR. These 
arguments are passed through either by OPPSL1 or OPEQEL. The 
Subroutine does the fol lowing: 


1, It identifies the cut generating element. 


2, The entire cut is never really appended to the tableau. Each 
column is updated by the pivot column in succession. The cut 
coefficient for each column is generated just before the column 
1s to be updated. 


3, Updating operations are straightforward. Adding or subtracting 


the pivot column multiplied by a scalar is all that is required to 
update a column. 
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Documentation of subroutine OPAUT I] 


Control enters OPAUT! either when optimality has been proved 
or when the user's iteration limit has been exceeded. It 
merely writes out the current solution to the problem. It 
also names the rows in order of the way they were input; 
$001,S002, etc. These are the names of the slack variables 
associated with each constraint. 
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ICETRAN LISTING OF SUBROUTINE OPAUT1 


24 


COMMON IAM(P 


Bue 
COMMON Y(P),E 


»B(P) ,NAMR(P) ,A(P),X(P) ,NAMC(P) , IROF(P) 
( ,P(P ie ICLF(P), KM(50) ,KBCD(9) ,KP(10) 
COMMON 2(15), T(6) KIT, IRW,LR,LK,LC,KPIV,KZZ,1S,IPZ,KZ 
COMMON NCAS (10 KL, IE,NCL,NRO sES(ONe IPIV, LEC 
DYNAMIC ARRAY IAM SNAMR,NAMC, IROF, ICLF 

DIMENSION IOPT(1) 

EQUIVALENCE (QQDUB(1),10PT(1)) 

WRITE(6,101) KPIV 

WRITE(6,103) IPZ 

EUECTS)Z4..23,24 

TRWO=IS-1 

DO 25 JAG=2,IRWO 

EEG O)=LES(5)+1 


B( 
P) 
(6 
) 


25 NAMR(JAG)=LES(5) 


23 
26 


nOZ 
103 
101 


GO TO 26 

PES (5 }=LES( 5) +1 

MAMR(IRW)=LES(5) 
11=(IFIND( IAM, 1)+QQICOM) /4 

DO 5 JxX=1,1RW 

JJZ=JX+11 

NRT=NAMR (JX) 

IZZ=I0PT(JJZ) 

WRITE(6,102) NRT,IZZ 

FORMAT(/,6X ,A4,'=',16) 

FORMAT(//,6X, ‘STATIONARY ITERATIONS’ ,2X,14) 
FORMAT(//,6X,'TOTAL ITERATIONS’, 2X,14) 
KZZ=-10 

RETURN 

END 
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Documentation of subroutine OPADRW 


Control enters OPADRW when the INSERT ROW command is encountered. 
It takes a row in the representation of the original problem 
Statement and updates 1t to the current tableau representation. 
This is don by multiplying the row by the tableau matrix that 
was previously put in by OPTAB as -I. When the new row has 

been expressed in the new tableau representation, control is 
returned to the OPCHOZ subroutine to see if the addition of this 
new row has changed the optimal solution to this problem. If 

it has, then the whole problem is reoptimized. 
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ICETRAN LISTING OF SUBROUTINE OPADR\ 


40 
35 


4] 


39 


35 
o7 


36 
38 


63 
62 


RO 


WwW 


SUBROUTINE OPADRW 

COMMON IAM(P),B(P),NAMR(P),A(P),X(P),NAMC(P),IROF(P) 
COMMON Y(P),E(P),P(P),ICLF(P) ,KM(50),KBCD(9) ,KP(10) 
COMMON Z(15),7T(6) KIT, IRW,LR,LK,LC,KPIV,KZZ,1S,IPZ,KZ 
COMMON NCAS(10) ,KL,TE,NCL,NRO,LES(10),IPIV,LEC 
DYNAMIC ARRAY IAM,NAMR,NAMC ,IROF, ICLF 

DIMENSION IOPT(1) 

EQUIVALENCE (QQDUB(1),IOPT(1)) 

IF(KZZ+10) 40,41,40 

WRITE(6,101) 

WRITE(6,104) 

CALL EXIT 

IcT=0 

IOP=NCL+2 

IEE=IRW+1] 

IF(LES(2)-LR) 35,39, 35 

LIT=] 

IROF(IEE)=1 

GO TO 38 

IF(LES(3)-LR)36, 37,36 

lIT=-] 

IROF(IEE)=0 


GO TO 38 

IROF(IEE)=-1 

IIT=] 

RELEASE IROF,LOW 

TZ0=INP-1 

DO 2 KX=1,I10P 

ICT=ICT+1 

READ(5,100) ITRM,NCCK,NRCK, ITEL 
DRUUCH]) 7.97 

WRITE(6,100) ITRM,NCCK,NRCK,ITEL 
IF( ITRM-LES(4) )3,613 

DO 62 KTR=1,1Z0 

IF (NAMC(KTR)-NCCK)62,63,62 
ICLF(KTR)=ITEL*IIT 

GO T0 2 

CONTINUE 

WRITE(6,100) ITRM,NCCK,NRCK,ITEL 
WRITE(6,105) NCCK 

GO TO 95 

CONTINUE 

WRITE(6, 102) 

GO TO 95 

IF(ICT-I0P)5,4,5 
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ICETRAN LISTING CONTINUED-OPADRW 


5 
4 
45 


44 
46 


47 
90 


100 
101 
102 
103 
104 
105 


WRITE(6,103) 

GO TO 95 

DO 50 JX=1,1Z0 

IF(JX-1)44,45,44 

TAM(JX, IEE) =ICLF (JX) 

GO JO 46 

TAM( JX, IEE )=0 

CALL ALOCAT (IAM, JX) 

11=(IFIND(IAM,JX)+QQICOM) /4 

133=11+IEE 

TOT=1 

DO 47 KX=IS,IE 

IOT=10T+1 

I] J=11+KX 

IOPT(133)=I0PT (133) -(1OPT(111)*ICLF(I0T) ) 

RELEASE IAM(JX),LOW 

CONTINUE 

RELEASE ICLF,LOW 

KPIV=0 

IPZ=0 

IRW=IEE 

LEC=-3 

CALL OPCHOZ 

RETURN 

FORMAT (3A4,2X,16) 

FORMAT (//,6X,*' INSERT ROW COMMAND OUT OF SEQUENCE. ') 
FORMAT (//,6X,'TOO MANY COEFF. IN INSERT ROW COMMAND. ') 
FORMAT(//,6X,'TOO FEW COEFF. IN INSERT ROW COMMAND. ') 
FORMAT(//,6X,'JOB FLUSHED. ') 

FORMAT(//,6X,'NO MATCH FOUND FOR-',A4,'-IN ABOVE CARD') 
END 
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